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ABSTRACT 


The  problem  of  time  difference  of  arrival  (TDOA)  is  important  in  underwater 
acoustics  for  both  passive  and  active  sonar.  Classical  approaches  to  this  problem  are 
based  on  generalized  cross-correlation  (GCC)  methods  implemented  in  the  frequency- 
domain.  After  appropriate  weighting  of  the  cross  spectral  data  in  the  frequency 
domain,  an  inverse  discrete  Fourier  transform  (IDFT)  is  performed  and  the  peak  of 
the  resulting  GCC  function  is  located  in  the  time  domain. 

This  thesis  shows  that  the  cross-spectrum  of  the  data  satisfies  an  appropriate 
signal  subspace  model;  therefore  the  IDFT  can  be  replaced  with  a  signal  subspace 
technique  such  as  MUSIC.  The  result  is  an  enhanced  abihty  to  locate  the  peak.  Fur¬ 
ther,  apphcation  of  methods  such  as  root-MUSIC  or  ESPRIT  produce  direct  numeri¬ 
cal  estimates  for  TDOA  without  the  need  to  search  for  a  peak.  Results  are  presented 
for  an  extensive  set  of  simulations  using  both  synthetic  signal  data  and  data  from 
a  ocean  acoustic  propagation  model  (MMPE).  Results  are  further  presented  for  an 
apphcation  of  the  new  method  to  target  localization  and  tracking.  In  all  cases  results 
are  compared  using  both  the  new  methods  and  the  classical  methods. 
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EXECUTIVE  SUMMARY 


In  the  passive  sonar  problem,  signals  received  at  two  or  more  sensors  or  hy¬ 
drophones  can  be  used  to  estimate  the  position  and  velocity  of  a  detected  acoustic 
source.  Passive  systems,  unlike  radar  or  active  sonar  systems,  cannot  control  the 
amount  of  transmitted  energy  to  be  reflected  from  the  source.  However,  the  covert¬ 
ness  of  passive  systems  can  be  advantageous,  both  in  military  and  biochemical  apph- 
cations.  In  practice,  the  number  of  receiving  sensors,  the  observation  time  and  the 
ratio  of  the  background  noise  to  the  source  signal  strength  after  propagation  loss, 
when  balanced  against  total  system  cost,  dictate  the  feasibility  of  passive  systems. 

In  the  ocean,  soxmd  usually  arrives  at  each  individual  omnidirectional  receiver 
through  more  than  one  path  (e.g.,  direct  and/or  surface,  bottom  reflected  paths).  In 
order  to  deal  with  this  problem  from  a  signal  processing  point  of  view,  it  is  useful  to 
decouple  two  separate  cases:  the  multipath  and  the  planar  problems.  For  multipath 
signals,  a  simple  model  of  the  received  signal  is  that  each  receiver  sees  a  signal  plus  an 
attenuated  and  delayed  signal  corrupted  by  additive  uncorrelated  noise,  coming  from 
different  directions.  For  the  planar  problem  (i.e.,  when  all  receivers  and  the  source 
are  in  the  same  plane),  it  is  usually  assumed  that  the  energy  arrives  at  each  receiver 
through  only  one  propagation  path  in  the  same  plane  with  all  receivers  and  source. 
This  means  that  the  time  delay  to  be  estimated  is  the  travel  time  of  the  acoustic 
wavefront  between  pairs  of  receivers,  so  that  the  source  position  and  velocity  can  be 
estimated. 

In  this  thesis  we  deal  with  the  time  difference  of  arrival  (TDOA)  problem. 
We  begin  by  examining  previous  research  for  estimating  time  difference  of  arrival 
known  as  generalized  cross-correlation  (GCC)  methods.  We  refer  to  these  as  “classical 
methods.”  With  these  methods,  data  from  the  two  channels  are  transformed  to 
the  frequency  domain  to  form  the  conjugate  product  XY*  (cross-spectrum).  After 
appropriate  weighting,  an  inverse  discrete  Fourier  transform  (IDFT)  is  performed  and 
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the  peak  of  the  resulting  generahzed  cross-correlation  function  is  located  in  the  time 
domain.  Our  main  goal  in  this  thesis  was  to  prove  that  for  this  kind  of  problem  it 
is  possible  to  use  so-called  signal  subspace  methods.  More  specifically,  we  showed 
that  the  aforementioned  cross-spectrum  of  the  data  satisfies  an  appropriate  signal 
subspace  model.  In  this  way  we  were  able  to  replace  the  IDFT  with  a  signal  subspace 
techmque  such  as  MUSIC,  Minimum  Norm,  PCLP,  or  similar  procedure.  In  our 
problem  we  apply  the  subspace  methods  in  the  firequency  domain  and  observe  their 
results  in  the  time  domain.  This  is  in  contrast  to  their  more  usual  apphcation  for  the 
estimation  of  signals  in  noise.  We  also  applied  the  subspace  methods  of  root-MUSIC 
and  ESPRIT,  which  produce  direct  numerical  estimates  for  TDOA,  without  the  need 
to  search  for  a  peak  in  the  “pseudo-correlation”  function.  Both  types  of  methods 
performed  efficiently. 

We  continued  with  a  thorough  series  of  simulations  using  synthetic  data  (MAT- 
LAB)  and  data  created  by  an  ocean  acoustic  propagation  model  (MMPE),  each  time 
testing  the  performance  of  all  methods.  The  results  showed  that  the  subspace  meth¬ 
ods  performed  as  well  or  better  than  the  classical  methods,  especially  with  low  SNR 
conditions  and  with  more  difficult  environmental  data.  Finally  we  applied  the  var¬ 
ious  TDOA  algorithms  (classical  and  subspace  methods)  to  the  problem  of  target 
localization  and  tracking.  All  methods  produced  successful  results. 
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I. 


INTRODUCTION 


The  problem  of  time-delay  estimation  (TDE)  is  important  in  the  field  of  imder- 
water  acoustics  for  both  passive  and  active  sonar.  In  this  problem  a  signal  is  emitted 
from  a  source  (e.g.,  a  submarine)  and  is  received  at  two  spatially  separated  sensors. 
If  si{t)  represents  the  original  undistorted  source  and  ni{t)  and  n^it)  represent  se¬ 
quences  of  imcorrelated,  additive  noise,  then  the  signals  received  at  two  spatially 
separated  sensors  may  be  modeled  as 

xi(t)  =  Si{t)  +  ni{t) 

X2{t)  =  asi{t  +  D)  +  n2{t)  (I.l) 

where  D  is  the  time  delay  between  the  two  sensors.  In  a  more  complicated  scenario, 
the  received  signals  are  attenuated  and  contain  multipath  reflections  of  the  original 
source  created  by  propagation  through  the  ocean  environment,  in  addition  to  the 
noise,  as  shown  in  Fig.  1.  In  the  rmknown  source  case  (generally  passive  sonar),  Xi(t) 
and  X2(t)  and  the  ordinary  cross  correlation  may  be  used  to  estimate  the  position 
and  velocity  of  a  moving  source.  For  the  known  source  case  (generally  active  sonar) 
ordinary  correlation  involves  using  the  original  source  and  only  one  received  signal 
to  estimate  the  time  it  takes  the  signal  to  travel  from  the  source  to  the  receiver,  i.e., 
ni(t)  =  0  in  the  model  given  above.  In  the  absence  of  propagation  distortion,  this  is 
matched  filtering,  which  is  the  optimum  linear  method  when  the  noise  is  white. 

Once  the  signal  has  been  detected  and  the  time  delay  D  has  been  estimated, 
the  time  delay  can  be  used  to  estimate  the  bearing  angle  B  shown  in  Fig.  2.  The 
bearing  estimate  is  given  by  the  approximate  rule 

B  ^  cos-\cD/L)  (1.2) 

where  c  is  the  speed  of  sound  in  water,  B  is  the  bearing  estimate  and  D  is  the  time 
delay  estimate.  It  can  be  shown  that  B  is  the  angle  that  the  hyperbohc  “line  of 
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Figure  1.  Model  of  direct  and  surface-reflected  sound  ray  paths  received  at  three 
sensors. 


Figure  2.  Planar  model  of  two  receivers  separated  by  a  distance  L. 

position  makes  with  the  axis  of  the  receivers;  hence  the  approximation  Eq.  (1.2)  is 
increasingly  accurate  as  the  range  to  the  acoustic  source  increases. 

The  critical  part  for  the  passive  bearing  estimation  problem,  or  in  general 
for  the  localization,  is  the  accurate  estimation  of  time  delay.  In  the  literature,  this 
problem  is  treated  using  terms  like  time  delay  estimation  (TDE),  time  delay  (TD), 
time  diflference  of  arrival  (TDOA),  group  delay,  time-of- arrival  diflerence  (TOAD), 
phase  delay  and  others.  For  our  purposes  it  is  assumed  that  all  these  terms  are 
identical  and  the  terms  TDE  or  TDOA  will  be  used  for  the  rest  of  our  discussion. 
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A.  PREVIOUS  RELATED  RESEARCH 


Two  conceptual  procedures  are  basically  known  according  to  [Ref.  1]. 

•  An  intuitive  approach  and  a  familiarity  with  detection  theory. 

•  A  rigorous  application  of  the  maximum  likelihood  (ML)  for  white  signals  in 
white  noise. 

In  both  conceptual  procedures  it  is  attempted  to  “advance”  the  delayed  received  sig¬ 
nal  by  a  hypothesized  amoimt  in  order  to  ahgn  it  with  the  other  received  signal.  In 
other  words  both  received  signals  Xi{t),  X2{t)  contain  the  same  transient  Si(t),  coming 
from  the  same  source,  “buried”  in  different  noise,  ni{t),  n2{t),  and  the  only  difference 
between  them,  is  the  time  of  arrival  (TOA).  After  the  above  hypothesis  either  they 
are  summed,  squared  and  averaged  as  shown  in  Fig.  3  or  multiphed  and  averaged 
as  shown  in  Fig.  4.  In  both  cases  the  hypothesized  delays  are  adjusted  in  order  to 


mautniz*  J(0) 

Figure  3.  Conceptual  delay,  sum,  square  and  integrate  configuration. 

maximize  the  configuration  output.  From  the  figures  it  is  observed,  that  both  config¬ 
uration  outputs  consist  of  “signal-cross-signal”  terms.  Further,  both  configurations 
are  ML  estimators  for  time  delay  imder  the  assumptions  that  the  signal  and  noises 
are  white  and  mutually  xmcorrelated  [Ref.  1].  When  the  signal  and  noise  spectral 
characteristics  are  nonwhite,  the  received  waveforms  must  be  prefiltered  with  partic¬ 
ular  equiphase  filters  (i.e.,  the  prefilters  must  have  the  same  phase  characteristics)  to 
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Figure  4.  Conceptual  cross  correlator  configuration. 


accentuate  the  frequency  bands  with  good  signal-to-noise  ratio  (SNR).  It  is  of  inter¬ 
est  to  note  that  a  signal  detector  can  be  realized  by  comparing  either  configuration 
output  to  a  threshold.  Moreover  in  terms  of  detection  (but  not  estimation)  in  the 
presence  of  noise,  the  system  in  Fig.  3  outperforms  the  system  in  Fig.  4  [Ref.  1]; 
however,  the  system  in  Fig.  3  requires  prior  knowledge  of  power  levels  in  order  to 
set  the  proper  threshold.  The  system  in  Fig.  4  has  a  zero  mean  output  in  the  signal 
absent,  noise  present  case. 

A  number  of  different  approaches  have  been  taken  so  that  the  conceptual 
systems  in  Figs.  3  and  4  can  be  achieved.  The  system  in  Fig.  3  can  be  implemented 
as  shown  in  Fig.  5  which  is  called  a  time  delay  beamformer.  It  is  presumed  that  s(t)  is 
a  sampled  version  of  the  original  broadband  time  signal  and  that  the  sampling  rate  is 
large  in  comparison  with  the  required  Nyquist  rate.  In  particular,  the  sampling  rate 
is  much  greater  than  twice  the  bandwidth.  The  system  configuration  shown  in  Fig. 
4,  including  prefilters  can  be  instrumented  by  a  generalized  cross-correlation  (GCC) 
function,  which  in  general  is  the  inverse  fast  Fourier  transform  (FFT)  of  the  product 
of  a  weighting  function  and  the  estimated  complex  cross-power  spectrum.  This  means 
that  data  from  the  two  sensors  are  transformed  to  the  frequency  domain  to  form  the 
conjugate  product  XY*.  After  appropriate  weighting,  an  inverse  DFT  is  performed 
and  the  peak  of  the  resulting  generalized  cross-correlation  function  is  located  in  the 
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Figure  5.  Time-domain  beamformer  implementation  of  one  conceptual  configuration 
for  three  hypothesized  delays. 

time  domain.  All  the  methods  that  follow  this  general  approach  will  be  referred  to  in 
this  thesis  as  classical  methods  and  quite  extensive  research  has  been  conducted  on 
these  methods  [Ref.  2,  3,  4,  1,  5]. 

More  recently,  techniques  based  on  higher  order  statistics  such  as  bicorrelation 
and  bispectra  have  been  introduced  to  the  problem  [Ref.  6,  7,  8,  9].  These  methods 
are  claimed  to  have  some  advantage  when  there  is  Gaussian  noise  present  which  is 
correlated  between  sensors.  In  particular,  since  higher  order  cumulants  of  Gaussian 
noise  are  theoretically  zero,  these  methods  are  practically  “blind”  to  Gaussian  noise. 

B.  THESIS  OBJECTIVES 

This  thesis  has  the  following  goals.  The  first  objective  is  to  research  and 
develop  an  algorithm  that  can  be  used  for  Time-Difference-Of- Arrival  (TDOA)  esti¬ 
mation  based  on  subspace  methods.  The  next  objective  is  to  evaluate  the  performance 
of  this  algorithm  by  comparing  it  with  the  respective  algortihm  based  on  generalized 
cross-correlation  (GCC)  methods  implemented  in  the  frequency  domain,  through  a 
thorough  set  of  simulations  using  synthetic  (MATLAB)  and  model-based  acoustic 
(MMPE)  data.  Finally  the  last  goal  in  this  thesis  is  to  develop  localization  and 
tracking  algorithms  implementing  the  above  methods,  and  to  evaluate  their  perfor¬ 
mance  in  target’s  extraction  data  (position-course-speed). 
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C.  THESIS  APPROACH 

In  this  thesis  we  suggest  an  alternative  to  the  classical  methods  based  on 
some  more  modern  approaches  to  spectral  analysis.  In  particular,  we  focus  on  meth¬ 
ods  based  on  the  signal  subspace  concept.  We  observe  that  the  data  product  XY* 
computed  in  the  classical  methods  satisfies  a  signal  subspace  model  and  therefore  a 
number  of  well-known  techniques  such  as  MUSIC  [Ref.  10],  ESPRIT  [Ref.  11]  and 
others  can  be  used  to  estimate  the  time  delay.  Mathematically,  the  procedure  can 
be  viewed  as  replacing  the  inverse  DFT  in  the  classical  methods  with  one  of  these 
other  techniques.  As  we  later  see  from  experimental  results  the  new  methods  retain 
the  accuracy  of  the  classical  methods,  but  give  a  much  “cleaner”  indication  of  the 
peak.  Further,  the  use  of  methods  such  as  root  MUSIC  or  ESPRIT  for  the  estimation 
provides  direct  nmnerical  estimates  for  the  time  delay  without  the  need  to  search  for 
a  maximum.  We  also  have  to  mention  that  we  have  obtained  similar  behavior  with 
the  “Maximum  Likefihood”  or  “Minimum  Variance”  method  proposed  by  J.  Capon 
[Ref.  12,  13].  Our  focus  in  this  thesis  however  is  on  the  subspace  methods. 

D.  THESIS  OUTLINE 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  II  provides  a 
short  introduction  of  the  concepts  of  the  classical  methods  and  also  describes  the 
special  characteristics  of  each  one  method,  that  is  going  to  be  used  later  in  the  sim¬ 
ulations  from  this  family  of  methods.  Chapter  III  discusses  the  general  idea  behind 
the  subspace  methods,  gives  a  brief  insight  of  all  the  methods  from  this  group,  and 
finally  provides  an  argumentative  presentation  of  the  way  that  these  methods  are  im¬ 
plemented  into  our  problem  (TDOA).  Chapter  IV  presents  the  whole  thesis  approach 
to  the  TDOA  problem,  using  various  synthetic  data  (MATLAB  software)  and  model 
based  acoustic  data  using  the  Monterey-Miami  Parabofic  Equation  (MMPE)  model 
[Ref.  14,  15]  for  more  reahstic  results,  under  different  conditions/environments  each 
time.  Chapter  V  discusses  the  simulation  results  for  both  kinds  of  data  (synthetic  and 
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model-based  acoustic)  for  the  TDOA  problem.  Chapter  VI  develops  a  localization 
algorithm  in  two  dimensions  and  discusses  the  implementation  of  the  above  methods 
in  this  locahzation  problem  in  order  to  examine,  how  they  ultimately  function  in  the 
desired  goal,  namely  the  “Target  Detection  and  Localization.”  Chapter  VII  further 
extends  the  research  to  “Target  Tracking  Problem.”  For  this,  one  has  to  estimate  not 
only  target  position,  but  also  target  course  and  speed,  using  Doppler  measurements 
provided  by  MMPE.  Chapter  VIII  presents  conclusions  and  suggestions  for  future 
work. 
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II.  CLASSICAL  METHODS 


A.  GENERAL  APPROACH 

We  have  seen  that  the  transmission  of  a  short  transient  signal  from  a  remote 
source  monitored  in  the  presence  of  noise  at  two  spatially  separated  sensors  can  be 
mathematically  modeled  as 

Xi{t)  =  Si{t)  +  ni{t)  , 

X2{t)  =  asi{t  +  D)  +  n2{t)  ,  (II.  1) 

where  Si{t),  ni{t),  and  n2(t)  are  real,  jointly  stationary  random  processes.  Signal 
Si{t)  is  assumed  to  be  uncorrelated  with  noise  ni{t)  and  n2(f). 

One  common  method  of  determining  the  time  delay  D  is  to  compute  the  cross 
correlation  function 

/  xi{t)x2it  -  T)dt ,  (II.2) 

where  T  represents  the  observation  interval.  In  order  to  improve  the  accuracy  of  the 
delay  estimate  D,  it  is  desirable  to  prefilter  Xi(t)  and  X2{t)  prior  to  integration  in  Eq. 
(II.2).  As  shown  in  Fig.  6,  Xi  may  be  filtered  through  Hi  to  yield  yi  for  i  =  1, 2.  The 
resultant  j/i  are  multiphed,  integrated,  and  squared  for  a  range  of  time  shifts,  r,  until 
the  peak  is  obtained.  The  time  shift  causing  the  peak  is  an  estimate  of  the  true  delay 
D.  When  the  filters  =  i?2(/)  =  IjV/,  the  estimate  D  is  simply  the  value  of  r 

at  which  the  cross-correlation  function  peaks. 

The  cross-correlation  between  xi{t)  and  X2{t)  is  related  to  the  cross-power 
spectral  density  function  by  the  Fourier  transform  relationship 

■Rx.»(r)  =  r  (11.3) 

»/— oo 

When  xi(t)  and  X2{t)  have  been  filtered  as  shown  in  Fig.  6,  the  cross-power  spectrum 
between  the  filter  outputs  is  given  by 

(II.4) 
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Figure  6.  Received  waveforms  filtered,  delayed,  multiplied,  and  integrated  for  a  vari¬ 
ety  of  delays  until  peak  output  is  obtained. 

where  *  denotes  the  complex  conjugate.  The  generalized  cross-correlation  between 
Xi(t)  and  X2{t)  is  defined  by 

W  =  /“  (f)e?^’^df,  (11.6) 

where 

Mf)  =  Mf)m(f)  (II.6) 

and  denotes  the  general  frequency  weighting. 

In  practice,  Gxixiif)  is  not  known  a  priori,  and  only  an  estimate  of 

^xiX2if)  can  be  obtained  from  finite  observations  of  Xi{t)  and  X2(t).  Consequently, 
the  integral 

(^)  =  Mf)Gx,x2  (IL7) 

is  evaluated  and  used  for  estimating  delay.  Indeed,  depending  on  the  particular  form 
of '0p(/)  and  the  a  priori  information,  it  may  also  be  necessary  to  estimate  tpg{f)  in 
Eqs.  (II.5)  and  (11.6).  For  example,  when  the  role  of  the  prefilters  is  to  accentuate 
the  signal  passed  to  the  correlator  at  those  frequencies  at  which  the  signal-to-noise 
(S/N)  ratio  is  highest,  then  ^g(f)  can  be  expected  to  be  a  function  of  signal  and  noise 
spectra  which  must  either  be  known  a,  priori  or  estimated. 
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B.  PROCESSOR  WEIGHTING  EVALUATION 


Before  continuing  to  describe  the  weighting  functions  that  are  used  in  the 
generalized  cross-correlation  (GCC)  in  practice,  it  is  informative  to  examine  first  the 
effect  of  processor  weightings  on  the  shape  of  under  ideal  conditions.  For 

models  of  the  form  of  (II.  1),  the  cross  correlation  of  xi(t)  and  X2{t)  is 

RxxX2{t)  =  OiRsisAr  -  D)  +  Rn^n^i'r)-  (II-8) 

The  Fourier  transform  of  (II.8)  gives  the  cross-power  spectrum 

OnM  =  (II.9) 

If  ni{t)  and  n2(t)  are  uncorrelated  {Gnm^if)  =  0),  the  cross-power  spectrum  between 
Xi{t)  and  X2(t)  is  a  scaled  signal  power  spectrum  times  a  complex  exponential.  Since 
multiphcation  in  one  domain  is  a  convolution  in  the  transformed  domain,  it  follows 
for  G„in2(/)  =  0  that 


Rxixiir)  =  aRs^sA^)  0  d(t  -  D).  (11.10) 

where  @  denotes  convolution. 

One  interpretation  of  (II.  10)  is  that  the  delta  function  has  been  spread  or 
“smeared”  by  the  Fourier  transform  of  the  signal  spectrum.  If  Si(t)  is  a  white  noise 
source,  then  its  Fomrier  transform  is  a  delta  function  and  no  spreading  taJses  place.  An 
important  property  of  autocorrelation  functions  is  that  |iZss(r)|  <  i?ss(0).  Equality 
will  hold  for  certain  r  for  periodic  functions.  However,  for  most  practical  apphcations, 
equality  does  not  hold  for  r  ^  0,  and  the  true  cross  correlation  (II.IO)  will  peak  at  D 
regardless  of  whether  or  not  it  is  spread  out.  The  spreading  simply  acts  to  broaden 
the  peak.  For  a  single  delay  this  may  not  be  a  serious  problem.  However,  when  the 
signal  has  multiple  delays,  the  true  cross  correlation  is  given  by 

R^.^A'r)  =  Rs,sAr)  0  -  A).  (11.11) 
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In  this  case,  the  convolution  with  Rsisi{t)  can  spread  one  delta  function  into  an¬ 
other,  thereby  making  it  impossible  to  distinguish  peaks  or  delay  times.  Under  ideal 
conditions  where  ^fjGxixiif)  —  should  be  chosen  to  ensure  a  large 

sharp  peak  in  Ry^y2{j^)  rather  than  a  broad  one  in  order  to  ensure  good  time-delay 
resolution.  However,  sharp  peaks  are  more  sensitive  to  errors  introduced  by  finite 
observation  time,  particularly  in  cases  of  low  signal-to-noise  S/N  ratio.  Thus,  as  with 
other  spectral  estimation  problems,  the  choice  of  '4’gif)  is  a  compromise  between  good 
resolution  and  stability. 

Having  all  the  above  in  consideration,  we  are  equipped  with  the  appropriate 
background  for  the  role  that  ^g{f)  is  going  to  play.  Thus,  let  us  examine  some 
generalizations  of  the  cross-correlation  function  individually. 

1.  Roth  Processor 

The  weighting  proposed  by  Roth  [Ref.  4],  namely  ^ 

yields 

Equation  (11.13)  estimates  the  impulse  response  of  the  optimmn  linear  (Wiener-Hopf) 
filter 

which  “best”  approximates  the  mapping  of  xi{t)  to  X2it).  If  ni(t)  5^  0,  as  is  generally 
the  case  for  (II.l),  then 


Gxixi(f)  —  Gsis^{f)  -1-  Gnimif), 


(11.15) 


and 


^yiyz 


(r)  =  -D)®  /” 


^*191  (/)  +  Gniuiif) 


(11.16) 


^the  subscript  R  is  used  here  to  distinguish  the  choice  of 
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Therefore,  except  when  Gnmiif)  equals  any  constant  (including  zero)  times  Gg^si{f)y 
the  delta  function  will  again  be  spread  out.  The  Roth  processor  has  the  desirable 
effect  of  suppressing  those  frequency  regions  where  Gnmiif)  is  large  and  GinjC/)  i® 
more  likely  to  be  in  error. 


2.  Smoothed  Coherence  Transform  (SCOT) 

Errors  in  GxiX2(/)  naay  be  due  to  frequency  bands  where  Gn2n2if)  is  large,  as 
well  as  bands  where  Gnmiif)  i®  iai’ge.  One  is  therefore  uncertain  whether  to  form 
’’Pnif)  =  w'  /  A  or  i^nif)  =  tt— TTT!  hence,  the  SCOT  [Ref.  3,  16]  selects 


*XiXi 


if) 


^X2X2{f) 

Mf)  ^ 


\/GxiXi  {f)G  X2X2  if) 

This  weighting  gives  the  SCOT  generalized  cross-correlation 

J  —OQ 


(11.17) 


(11.18) 


where  the  coherence  estimate  is  given  by 


For  Hiif)  = 


\jGxixiif) 


^XlX2if)  — 
and  H2if)  = 


GxiX2if) 


iJGxiXi  (/)G,.„(/)' 

1 


(U.19) 


yj  Gx2X2if) 


,  the  SCOT  can  be  interpreted  (see 


Fig.  6)  as  prewhitening  filters  followed  by  a  cross  correlation.  When  Gxmiif)  = 
Gx2X2if)i  the  SCOT  is  equivalent  to  the  Roth  processor.  If  ni  0  and  722  7^  0,  the 
SCOT  exhibits  the  same  spreading  as  the  Roth  processor.  This  broadening  persists 
because  of  an  apparent  inabihty  to  adequately  prewhiten  the  cross  power  spectnnn. 


3.  Phase  Transform  (PHAT) 

To  avoid  the  spreading  evident  above,  the  PHAT  uses  the  weighting  [Ref.  3] 


which  yields 


i’pif)  = 


1 

\Gxix2if)y 


GxiX2if) 

\GxiX2if)\ 


(11.20) 


(11.21) 
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For  the  model  (II.l)  with  imcorrelated  noise  (i.e.,  Gmmif)  =  0), 


l^?x,x,(/)|  =  aGs^sAf)- 


(IL22) 


Ideally,  when  =  G^,^^{f), 

^xixaC/)  _  j2nfD 

has  unit  magnitude  and 

R^^L(.r)  =  5(t  -  D). 


(11.23) 


(11.24) 


The  PHAT  was  developed  purely  as  an  ad  hoc  technique.  Notice  that  for 
models  of  the  form  of  (II.l)  with  uncorrelated  noises,  the  PHAT  (11.21)  ideally  does 
not  suffer  from  the  spreading  that  other  processors  do.  In  practice,  however,  when 
Gxixiif)  ^  Gx-^x^if)-!  then  d{f)  ^  2'KfD  and  the  estimate  of  Ryfy^i'^)  will  not  be  a 
delta  function.  Another  apparent  defect  of  the  PHAT  is  that  it  weights  Gx^x^if)  as 
the  inverse  of  Gs^si  (/)•  Thus  errors  are  accentuated  where  signal  power  is  smallest.  In 
particular,  if  Gx^x^if)  =  0  in  some  frequency  band,  then  the  phase  ${f)  is  undefined 
in  that  band  and  the  estimate  of  the  phase  is  erratic,  being  imiformly  distributed  in 
the  interval  [-7r,7r]  radians.  For  models  of  the  form  of  (III),  this  behavior  suggests 
that  'ipp{f)  should  be  additionally  weighted  to  compensate  for  the  presence  or  absence 
of  signal  power.  The  SCOT  is  one  method  of  doing  this. 


4.  The  Eckart  Filter 

The  Eckart  filter  derives  its  name  from  work  in  this  area  done  in  [Ref.  17]. 
Although  it  is  not  used  in  the  experimental  part  of  the  thesis,  derivations  in  [Ref. 
18,  19,  20]  and  [Ref.  21]  are  outlined  here  briefly  for  completeness.  The  Eckart  filter 
maximizes  the  deflection  criterion,  i.e.,  the  ratio  of  the  change  in  mean  correlator 
output  due  to  signal  present  to  the  standard  deviation  of  correlator  output  due  to 
noise  alone.  For  long  averaging  time  T,  the  deflection  has  been  shown  [Ref.  19]  to  be 

roc  \Hl{f)m2(fWGn,nAf)Gn,nMdr  ^ 
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where  L  is  a  constant  proportional  to  T,  and  GsiS2(f)  is  the  cross-power  spectrum 
between  si{t)  and  s^it).  For  the  model  (II.  1),  Gais^if)  =  a:(j3isi(/) exp(j27r/D). 
Application  of  Schwartz’s  inequality  to  (11.25)  indicates  that 

HiitWif)  =  (11.26) 


maximizes  <P  where 


aG 


SiSi 


(/)G  (/)■ 


(11.27) 


Notice  that  the  weighting  (11.27),  referred  to  as  the  Eckart  filter,  possesses  some  of 
the  qualities  of  the  SCOT.  In  particular,  both  weightings  act  to  suppress  frequency 
bands  of  high  noise.  Also  note  that  the  Eckart  filter,  unlike  the  PHAT,  provides  us 
with  zero  weight  when  Gg^si  =  0.  In  practice,  the  Eckart  filter  requires  knowledge 
or  estimation  of  the  signal  and  noise  spectra.  For  (II.l),  when  o:  =  1  this  can  be 
accomplished  by  letting 


Mf)  =  -  |4.»(/)I1  •  [G«,.(/)  - 


(11.28) 


Since  in  our  problem  it  was  not  possible  to  have  knowledge  or  estimation  of  the 
signal  and  noise  spectra,  it  was  not  possible  to  use  this  weighting  processor  for  the 
experimental  part  of  the  thesis.  All  the  above  processors  are  summarized  in  Table  I 
and  can  be  justified  on  the  basis  of  reasonable  performance  criteria,  whether  heuristic 
or  mathematical. 


Processor  Nzune 

Weight 

M)  = 

Cross  Correlation 

1 

Roth  Impulse  Response 

llGxiXi{f) 

SCOT 

1/  y  GxiXi  {f)Gx2X2if) 

PHAT 

l/\GxiX2{f)\ 

Eckart 

^3iSi/[Gnini  (/)G^7^27^2  (/)] 

Table  I.  Candidate  Processors 
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III.  SUBSPACE  METHODS 


A.  CONCEPT  AND  PRINCIPLES 


An  entire  class  of  spectral  estimates  is  based  on  the  concept  of  signal  and  noise 
subspaces  associated  with  the  correlation  matrix  for  a  random  process.  Subspace 
methods  apply  primarily  to  locating  discrete  components  (hnes)  of  the  spectrum. 
Pisarenko’s  harmonic  decomposition  was  the  first  of  these  methods  and  motivated 
other  improved  methods,  such  as  MUSIC,  that  followed.  These  methods  are  all 
based  on  the  fact  that  when  the  data  consists  of  complex  exponentials  in  noise,  the 
frequency  vector  w,  defined  by 


1 


(111.1) 


is  an  eigenvector  of  the  correlation  matrix.  Its  projection  onto  an  orthogonal  subspace 
complementary  to  the  subspace  defined  by  all  of  the  signal  vectors  produces  a  null 
which  can  be  exploited  to  estimate  the  frequency.  The  methods  can  be  formulated 
either  as  a  search  for  peaks  in  a  function  (called  a  pseudospectrum)  or  as  a  polynomial 
root-finding  problem.  In  this  section  we  give  a  short  presentation  of  the  general 
principles  [Ref.  22]  in  these  methods. 

Consider  M  independent  signals  in  noise,  the  observation  vector  x,  the  noise 
vector  *7,  and  the  signals  vectors  Si,  defined  by  N  consecutive  samples  of  the  random 
process,  and  where  it  is  assumed  that  M  <  N.  Analytically,  this  means  that  the 
observation  signals  x[n]  (transient  Sz  noise)  are  equal  with 

M 

=  £  SiN  +  vln]  ,  (III-2) 

f=l 

where  the  transients  Si[n]  are  given  by 


Si[n]  =  . 


(III.3) 
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Specifically, 


x[0] 

x[l] 

.  >  (III.4) 

x[N  —  1] 

77[0] 

,  ^[1] 

^  =  I  .  >  (IIL5) 

r?[iV  -  1]  _ 

1 

(IIL6) 

x  =  53AiSi  +  r;,  (III7) 

i=l 

where  Si  is  defined  by  (III.6),  and  Ai  is  the  complex  amphtude  of  the  signal, 

A  =  (IIL8) 


This  is  the  general  signal  model  used  in  all  of  the  subspace  methods.  If  the  noise  is 
white,  the  correlation  matrix  is 


M 

Rx  =  PiSiS^  ^  +  a|l ,  (IIL9) 

i=l 

where  Pi  is  defined  by 

Pi  =  E{AiA*}  =  E[\Ai\^]  .  (III.IO) 

The  last  two  equations  can  also  be  written  with  more  compact  matrix  notation  as 

Ai 
A2 

+  »7  (III.ll) 

Am 
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and 


=  SPoS*'^  +  all , 


(III.  12) 


where 


and 


S  = 


Si  S2  •  •  •  Sm 


Po 


Pi  0  •  •  •  0 
0  P2  •••  0 


(III.13) 


(III.  14) 


[0  0  ■■■  Pm  \ 

It  can  be  shown  that  the  correlation  matrix  Rx  has  M  eigenvectors  lying  in  the  signal 
subspace,  which  is  spaimed  by  the  Si,  and  N—M  eigenvectors  lying  in  the  orthogonal, 
complementary  noise  subspace.  All  TV  —  M  noise  subspace  eigenvectors  correspond 
to  eigenvalues  Xi  =  al  while  all  of  the  signal  subspace  eigenvectors  correspond  to 
eigenvalues  Xi>  al  . 

When  the  noise  is  not  white,  but  is  still  uncorrelated  with  the  signals,  the 
foregoing  development  leads  to  the  correlation  matrix 


=  SPoS*'^  +  al^r, 


(III.  15) 


instead  of  (III.  12).  Here  S,,  is  a  normahzed  covariance  matrix  that  represents  the 
covariance  structure  of  the  noise  vector.  In  this  case  the  whitening  transformation 

y  =  (III.16) 


leads  to  the  correlation  matrix 


Ry  =  =  TPoT*"^  +  all  ,  (III.17) 


where 

T  =  (III.  18) 
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is  the  matrix  with  columns  that  are  the  transformed  signal  vectors 

tk  =  S-^/^Sk;  k  =  .  (111.19) 

In  the  transformed  vector  space,  the  eigenvectors  corresponding  to  the  M 
largest  eigenvalues  span  the  signal  subspace;  those  corresponding  to  the  smallest 
eigenvalues  (equal  to  a^)  span  the  noise  subspace.  The  eigenvectors  and  eigenvalues 
satisfy  the  equation 

Hyelc  =  =  Ake^  (III.20) 

which  can  be  written  as  the  generaUzed  eigenvalue  problem 


I^®k  —  AkS,jek  , 


(III.21) 


where 

ek  =  .  (III.22) 

The  basis  vectors  that  span  the  signal  and  noise  subspaces  in  the  original  coordinate 
system  are 

bk  =  (Sj/^'ek)  (III.23) 

or 


bk  =  S„ek  .  (III.24) 

Neither  the  eigenvectors  ek  nor  the  basis  vectors  bk  are  orthonormal  in  the  usual 
sense.  The  eigenvalues  satisfy  the  same  conditions,  however.  Those  corresponding 
to  the  noise  subspace  are  the  N  —  M  smallest  eigenvalues  (equal  to  a^)  and  those 
corresponding  to  signal  subspace  are  the  M  largest  eigenvalues  (all  greater  than  al). 

Let  us  now  return  to  the  case  of  white  noise  and  the  correlation  matrix  (III.  12). 
It  is  convenient  for  our  later  discussion  to  define  the  matrices  of  eigenvectors 


Esig  = 


Cl  62  •  •  •  eM 


(III.25) 
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and 


Enoise  — 


GM+I  eM+2 


Cn 


(III.26) 


and  also  the  two  matrices  of  eigenvalues 

Ai  0 


Asig  — 


0  As 


0  0 


0 

0 

Ajw 


(III.27) 


Am+1 

0 

...  0 

'  <^1  0 

...  0 

0 

Am+2 

...  0 

0  al 

...  0 

0 

0 

0  0 

•  •  •  (^1  _ 

and 


Anoise  — 


The  complete  matrices  of  eigenvectors  and  eigenvalues  are  thus  given  by 

E  —  [Egig  Enoise] 

and 

A= 

0  1 

Now  observe  that  it  is  possible  to  write  Rx  as 

Rx  =  EAE*""  =  EsigAsigE^  +  EnoiseAnoiseE^^i^, 

and  to  write  R^^  as 

R;*  =  EA-'E-’’  =  E.,.Ar,jE:5  +  E„<„„A;^^ES^  . 


(III.28) 


0 


(111.29) 


(III.30) 


(IIL31) 


(111.32) 


Further,  since  the  columns  of  Egig  are  orthonormal  and  define  the  signal  subspace, 
this  matrix  can  be  used  to  form  a  projection  matrix  for  the  signal  subspace, 

,-i . 


_  TTl  /'C'^TtTI  \  T7‘^T  TTl 

sig  -^sig  ^-^sig-^sig  j  “^sig  -^sig-l^sig  > 


(III.33) 
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where  the  last  step  follows  because  E^jEsig  =  Imxm-  Likewise,  Enoise  can  be  used 
to  form  a  projection  matrix  for  the  noise  subspace 

Pnoise  =  E„oiseE;^i3^  =  I  -  ,  (III.34) 

where  the  last  equahty  follows  because  the  subspaces  are  orthogonal  complements  of 
each  other. 


B.  TYPES  OF  SUBSPACE  METHODS 
1.  MUSIC 


The  MUSIC  (for  Multiple  Signal  Classification)  method  forms  a  correlation 
matrix  of  some  size  N  >  M  +  1  from  which  we  try  to  find  the  eigenvalues  and 
eigenvectors.  Prom  the  previous  presentation,  the  eigenvalues  and  the  corresponding 
eigenvectors  are  divided  into  two  categories;  the  ones  that  “belong”  to  signal  subspace 
and  the  others  to  noise  subspace  according  to  the  estimated  eigenvalues.  If  the  number 
of  signals  is  not  known,  it  can  be  estimated  by  looking  at  the  smallest  eigenvalues 
and  finding  the  set  that  are  approximately  equal.  This  number  is  equal  to  N  -  M. 
The  MUSIC  method  involves  projection  of  the  signal  onto  the  entire  noise  subspace. 


For  the  MUSIC  method  we  have  the  following  approaches: 


1.  Use  the  squared  magmtude  of  the  projection  of  w  onto  the  noise  subspace. 
Since  each  of  the  signals  is  orthogonal  to  the  noise  subspace,  the  quantity 
w  PnoiseW  =  ■w*^EnoiseE*QjggW  goes  to  zero  for  the  values  of  the  frequency 
where  w  =  Si.  The  MUSIC  pseudospectrum  is  defined  as 


1 


1 


W*’'PnoiseW 

and  therefore  exhibits  sharp  peaks  at  the  signal  frequencies  where  w  =  Si. 


w*'rE„oiseE;iJi^w 


(III.35) 


2.  Use  an  alternative  root-finding  variation  of  the  method  called  “roof  MUSIC”. 
In  this  approach  an  eigenfilter  Ei{z)  is  defined  as 

Ei{z)  =  ei[0]  +  ei[l]z-^  +  •■■  +  e^AT  -  ,  (III.36) 

where  the  ei[n]  are  components  of  the  eigenvector  ei.  In  this  way  the  MUSIC 
pseudospectrmn  can  be  expressed  as 


(111.37) 
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Since  the  denominator  goes  to  zero  at  z  =  =  1,2,- 

nator  pol3aiomial 

N 


Pm\;(^)=  E  b,(z)e;{i/z') 


i=M+l 


•  • ,  M),  the  denomi- 
(III.38) 


has  M  roots  lying  on  the  unit  circle.  These  M  roots  (which  are,  in  fact,  double 
roots)  correspond  to  the  signal  frequencies. 


A  modification  to  the  MUSIC  method  was  proposed  by  Johnson  and  DeGraff 
[Ref.  23].  Prom  the  previous  discussion,  the  MUSIC  pseudospectriun  can  be 
written  as 


w' 


■»T 


(E£M+ieier’')w 


(III.39) 


Since  the  eigenvalues  for  all  the  noise  subspace  eigenvectors  are  equal  to  cr^,  a 
pseudospectrum  for  MUSIC  that  differs  firom  (III.35)  only  by  a  constant  can 
be  defined  as 


P MuW'^)  — 


1 


w 


*T 


(sr.: 


M+iejer-l  w  w 


♦T 


(SHm+I 


w 


(III.40) 


where  the  last  equahty  follows  because  \i  =  al,i  —  M  1,  -  •  •  .  Equation 

(III.40)  gives  the  pseudospectrum  of  the  “modified  MUSIC  method^’  by  John¬ 
son  and  DeGraff  which  differs  from  the  first  approach  in  that,  in  practice,  the 
estimated  eigenvalues  are  not  all  exactly  equal. 


2.  Minimum-Norm  Procedure 

In  the  Minimum-norm  procedure  [Ref.  24,  25]  we  find  a  single  appropriately 
chosen  vector  d  in  the  noise  subspace  and  define  the  pseudospectrum  in  terms  of  this 
vector  as 

=  |,^*Td|2  =  w*Tdd*Tw  ■ 

The  vector  which  lies  in  the  subspace  is  chosen  so  that  the  squared  magnitude  ||d|p 
is  minimized  subject  to  the  constraint  that  its  first  component  is  equal  to  1.  The 
resulting  vector  is  given  by 


d  = 


-E. 


C*^C 


(111.42) 


23 


where  Enoise  is  the  matrix  of  noise  subspace  eigenvectors  and  c*^  is  the  top  row  from 
the  partition  of  Enoise  such  that 


Enoise  — 


,*T 


El 


(III.43) 


If  the  components  of  d  are  denoted  by  d[0],  d[l],  1],  with  d[0]  =  1,  then  the 

minimum-norm  pseudospectrum  can  be  equivalently  expressed  as 


=  \D^  ’ 

where  D{z)  is  the  polynomial 

D{z)  =  ^  d[k]z-'^  .  (III.45) 

k=0 

Thus  the  frequencies  can  also  be  found  as  the  roots  of  D(z)  on  the  unit  circle. 

3.  Principal  Components  Linear  Prediction 

Another  technique  involving  simple  modifications  to  some  of  the  previously 
presented  methods  suggests  itself.  Whenever  there  is  added  noise,  let  the  estimated 
correlation  matrix  be  represented  in  terms  of  only  its  eigenvectors  and  eigenvalues 
pertaining  to  the  signal  subspace. 

M 

=  (111.46) 

i=l 

This  is  sometimes  referred  to  as  the  principal  components  approximation  of  the  core¬ 
lation  matrix.  This  procedure,  developed  by  Tufts  and  Kumaresan  [Ref.  26,  27], 
amounts  to  ehminating  some  of  the  noise  and  increasing  the  overall  signal-to-noise 
ratio.  For  model-based  methods,  a  model  derived  from  the  reduced  rank  correlation 
matrix  alone  can  then  be  generated  and  used  to  estimate  the  spectrum. 

Tufts  and  Kumaresan  tried  to  eUminate  the  terms  in  the  noise  subspace  using 
the  pseudoinverse.  This  task  could  have  even  more  successful  results  if  a  high-order 
correlation  matrix  was  implemented.  If  the  value  of  the  prediction  order  variable  P 


24 


was  approximately  in  the  same  range  with  the  available  data  samples  (Ng),  then  the 
rank  of  the  estimated  correlation  matrix  is  automatically  reduced.  In  order  to  have  an 
exact  equality  between  the  rank  of  the  estimated  correlation  matrix  and  the  number 
of  signals,  then  the  following  relationship  should  hold:  P  =  Ng  —  M/2,  where  P  and 
Ng  are  defined  above  and  M  is  the  rank  of  the  pseudoinverse.  This  specific  case  is 
called  Kumaresan-Prony  case.  It  has  the  advantage  of  less  computational  demands 
since  the  SVD  of  the  square  matrix  is  the  same  as  its  eigenvalue  decomposition,  and 
there  is  no  need  for  both  estimations  to  be  conducted.  For  minimizing  the  variance 
of  the  firequency  estimates,  however,  a  smaller  value  of  P  »  jNg  is  suggested  [Ref. 
26,  27].  The  filter  coeflicient  vector  is  thus  computed  from 


M  /  '*T„\ 


(III.47) 


where  and  r  are  the  appropriate  partitions  of  the  correlation  matrix  Rx,  and  the 
e[,  are  the  eigenvectors  and  eigenvalues  of  R^-  The  pseudospectrum  is  given  by 


the  final  expression 


where 


PpcLp{^‘^) 


|w*'^aP 


(III.48) 


(III.49) 


If  the  number  of  signals  M  is  not  known,  it  has  to  be  estimated  from  the  eigenvalues. 
The  signal  frequencies  are  then  determined  by  either  finding  the  peaks  of  (III.48)  or 
by  finding  the  roots  of  A{z). 


4.  ESPRIT 

The  ESPRIT  (Estimation  of  Signal  Parameters  via  Rotational  Invariance 
Techniques)  [Ref.  28,  29]  method  takes  a  somewhat  diflFerent  approach  to  frequency 
spectrum  estimation  by  exploiting  a  certain  invariance  principle  that  exists  for  the 
subspaces  of  two  overlapping  data  sets.  Let  us  consider  a  data  set  of  iV  +  1  samples 
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x[0],  a;[l],  •  •  • ,  2:[A/'^]  and  define  the  two  vectors 


a:[0] 

■  x[l]  ■ 

X  = 

^[1] 

and  x  = 

x[2] 

r 

1 

* 

1 _ 

(IIL50) 


Let  B  and  B  denote  two  matrices  whose  columns  are  basis  vectors  for  the 
vector  spaces  associated  with  x  and  x'.  It  can  be  shown  [Ref.  22]  that  B  and  B'  are 
related  by  a  set  of  hnear  equations  of  the  form 


=  B\  (III.51) 

where  ^  is  an  M  x  M  square  matrix  and  the  eigenvalues  of  ^  are  of  the  form  Ai  =  , 

where  oji  are  the  desired  frequencies. 

In  the  next  few  lines  there  will  be  a  short  presentation  of  the  steps  in  the  ES¬ 
PRIT  algorithm  (TLS  version)  m  order  to  have  an  understanding  of  its  functionality. 
We  start  by  defining  the  N  -f-  1-dimensional  random  vector  x  pertaining  to  W  -H  1 
consecutive  data  samples  a;[0],a:[l],  •  •  •  ,a:[W]  and  estimating  the  correlation  matrix 
R®  from  the  data.  Then  we  compute  the  eigenvectors  and  eigenvalues  of  R.^  :  ^ 

^fc  =  Afc2Aj€k  k  =  ,N  +  1  .  (III.52) 


If  necessary,  we  estimate  the  number  of  signals  M  from  the  eigenvalues.  Afterwards 
a  basis  spanning  the  signal  subspace  is  generated  and  is  partitioned  as 


(III.53) 


Since  B  and  B'  will  not  satisfy  III.51  exactly,  a  Total  Least  Squares  (TLS)  approach 
to  computing  is  used.  The  details  of  TLS  are  beyond  the  scope  of  this  presentation. 

^If  the  noise  is  not  white  III.52  can  be  replaced  by  a  generalized  eigenvalue  problem. 
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However,  a  general  outline  would  be  to  first  compute  the  matrix  V  of  right  singular 
vectors  of 

B  B' 

and  then  partition  V  into  four  M  x  M  submatrices 


(111.54) 


V  = 


Vn  Vi2 

V21  V22 


The  desired  estimate  for  ^  is  then  given  by 


(III.55) 


^tls  =  -V12V22"  . 


(III.56) 


Finally,  the  desired  frequencies  are  computed  as  the  angle  of  the  complex  eigenvalues 
of  ^tls  as  shown  below 


Uk  =  lXk]  =  (III.57) 

C.  APPLICATION  TO  THE  TDOA  PROBLEM 

A  model  for  the  TDOA  problem  as  we  have  already  seen  is  as  follows.  We 
assume  that  a  short  transient  signal  d[n]  is  emitted  from  the  source.  We  further 
assume  the  signal  received  at  each  sensor  is  subject  to  amplitude  attenuation  and 
additive  noise  so  that  the  signals  x  and  y  can  be  written  in  discrete  time  as 

x[n]  =  d[n]  +  w[ri\ 

y  [n]  =  Ad[n  —  L]  +  u[n] ,  (III.58) 

where  A  is  the  relative  amphtude,  L  is  the  time  delay  to  be  estimated,  and  u  and  w 
are  the  additive  noise  terms.  Note  that  we  assume  the  transient  signals,  d,  received 
at  each  location  axe  identical  but  arrive  at  different  times.  This  assmnption  will  be 
relaxed  in  the  simulations  presented  later  [Ref.  30] .  Let  us  define  Rffcj  as  the  product 

R[k]  =  X[A:]y*(fc],  (III.59) 
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where  and  are  the  DFT’s  of  the  sequences  x[n]  and  y\p\  computed  at  a 

size  equal  to  at  least  twice  the  length  of  the  data.  Thus  the  inverse  transform  of  i?[A:] 

represents  an  estimate  of  cross-correlation  in  the  time  domain.  Now  let  ujo  represent 
the  frequency  resolution  of  the  DFT.  The  DFT  sequences  X[A:]  and  y[A:]  are  then 
given  by 

X[k]  =  D[k]  +  W[k] 

Y[k]  =  AD[k]e-^‘^°’^^  +  U[k]  (111.60) 

and  from  (III.59)  and  (III.60) 

R[k]  =  -h  {D[k]U*[k]  +  AD*[k]W[k]€^‘^^^^  +  W[k]U*[k]}.  (III.61) 

Let  us  also  assume  that  the  signal  d[n]  is  broadband,  and  that  the  source  spectrum 
is  flat  and  therefore  the  magnitude  |D[fc]|  is  approximately  constant.  Further,  if 
the  observed  time  sequences  {x  and  y)  are  sufficiently  long,  the  points  in  the  DFT 
sequence  for  the  term  enclosed  in  brackets  will  be  approximately  uncorrelated.  Thus 
the  sequence  i?[fc]  satisfies  the  basic  model  for  signal  subspace  analysis  [Ref.  22],  but 
with  time  and  frequency  interchanged. 

For  the  following  discussion,  we  shall  consider  both  related  approaches  to 
estimating  TDOA.  In  the  first  approach  a  steering  vector  of  the  form 

is  projected  onto  the  noise  subspace  and  the  reciprocal  of  the  result  is  plotted  as 
the  lag  parameter  I  is  varied.  (Here  N  is  the  dimension  of  the  observation  used  in 
the  chosen  subspace  method.)  The  resulting  function,  which  we  refer  to  as  a  delay 
indicator  function  (DIF),  peaks  at  the  desired  estimate  I  =  L.  Although  the  DIF 
resembles  a  generalized  cross-correlation  function,  it  is  strictly  speaking  not  a  GCC. 
It  is  analogous  to  the  pseudospectrum  used  in  the  spectral  estimation  problem.  In 
the  second  approach,  direct  numerical  estimates  are  computed  and  no  DIF  is  formed. 
Root  MUSIC  and  ESPRIT  follow  this  approach. 
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IV.  THESIS  APPROACH  TO  THE 

PROBLEM 


In  order  to  verify  aU  the  theoretical  results  described  in  Chapters  2  and  3,  var¬ 
ious  simulations  were  conducted.  AU  the  simulations  are  divided  into  two  categories 
depending  on  the  kind  of  the  data  that  were  used.  The  first  category  used  synthetic 
data  generated  from  equations  in  MATLAB  while  the  second  category  used  data  gen¬ 
erated  by  the  Monterey-Miami  Parabohc  Equation  (MMPE)  acoustic  propagation 
model  [Ref.  14,  15]. 

AU  the  simulations  using  the  synthetic  data  assumed  the  same  shape  and  form 
and  length  of  transient  at  each  receiver  but  for  each  one  we  added  different  noise  with 
variations  according  to  kind  (white-colored),  characteristics  (variance),  and  amount  of 
noise  (SNR).  White  or  colored  Gaussian  noise  was  generated  in  MATLAB  and  added 
to  the  signals  to  achieve  a  specified  signal- to-noise  ratio  (SNR).  The  SNR  specificaUy 
was  defined  by  the  formula 

±rpN-l  2r„l 

SNR  =  — Li  ^  (IV.l) 

where  N  is  the  number  of  the  samples  of  the  signal,  s[n],  and  is  the  variance  of 
the  additive  noise  (white  or  colored).  For  the  MMPE  data,  the  received  transients 
were  different,  in  general,  and  noise  was  added  with  variations  as  before.  In  this  way 
we  kept  aU  the  basic  assumptions  that  we  made  in  Chapter  1  for  the  TDOA  problem. 


A.  SYNTHETIC  DATA 

For  the  synthetic  data  five  different  kinds  of  transients  were  used: 

1.  exponential  transient  of  the  form 

gu\  ^  [  exp(-dt),  0<t<T 
^  '  I  0,  otherwise  ’ 


(IV.2) 


where  d  is  the  decay  factor  and  T  is  the  time  duration  (length)  of  the  expo¬ 
nential  transient; 
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2.  sinusoidal  transient  of  the  form 


s{t) 


sin(27r/t),  0  <t  <T 
0,  otherwise  ’ 


(IV.3) 


where  /  is  the  frequency  and  T  is  the  time  duration  (length)  of  the  sinusoidal 
transient; 

3.  damped  sinusoidal  transient  which  is  a  combination  of  the  first  two  transients, 


s{t)  = 


exp(-dt)  •  sin(27r/t),  0  <  t  <  T 
0,  otherwise 


(IV.4) 


4.  linear  swept-frequency  cosine  generator  {chirp{T,  fo,Ti,  /i))  transient, 


s{t)  =  (  0<t<T 

[  0,  otherwise 


(IV.5) 


This  function  creates  samples  of  a  hnear  swept-frequency  cosine  signal  at  the 
time  instances  defined  in  array  T,  /o  is  the  instantaneous  frequency  at  time  0, 
and  /i  is  the  instantaneous  frequency  at  time  Ti.  Both  fo  and  fi  are  in  Hertz. 
In  our  case  the  instantaneous  frequency  sweep  fi{t)  given  by 


where 


fo  d"  pt,  0  <  t  <  T" 

0,  otherwise  ’ 


/?  = 


/i  ~  fo 
Ti  ’ 


(IV.6) 


(IV.7) 


5. 


damped  chirp  transient  which  is  the  combination  of  the  first  and  fourth  tran¬ 
sients. 


s(t)  = 


exp(-dt)  •  0<t<T 

0,  otherwise 


(IV.8) 


In  this  way  we  were  able  to  compare  the  performance  of  the  classical  and 
subspace  methods  not  only  by  varying  the  characteristics  of  the -noise  but  also  by 
changing  the  characteristics  of  the  original  transients.  This  was  done  in  order  to  have 
a  more  complete  evaluation  of  the  above  methods. 
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B.  MODEL  BASED  ACOUSTIC  DATA 


Although  they  provide  a  way  to  test  algorithms  imder  ideal  conditions,  we 
understand  the  simulated  data  in  MATLAB  does  not  approach  the  real  conditions 
that  occur  in  the  ocean.  For  this  reason  we  used  data  from  the  MMPE  model  to 
test  the  TDOA  methods  under  more  realistic  conditions.  This  allowed  us  to  test 
them  in  more  complex  environments,  where  additional  factors  like  area  geography, 
water  column  and  bottom  characteristics,  sound  speed  profile,  etc.,  can  change  and 
influence  the  result  of  the  sound  wave  propagation. 

1.  Monterey-Mieimi  Parabolic  Equation  (MMPE)  Model 
and  Split-Step  Fourier  (PE/SSF)  Algorithm 

Before  presenting  data  generated  by  the  MMPE  model,  it  is  informative  to 
describe  what  the  model  is  what  the  model  is  and  what  it  is  capable  of.  We  begin  by 
defining  the  time- harmonic  acoustic  field  in  cylindrical  coordinates  (r,  z,  (f)) 

P{r,  z,  0,  u;t)  =  p(r,  z,  (IV.9) 


The  result  of  the  combination  of  Eq.  (IV.Q)  with  the  wave  equation,  also  in  cylindrical 
coordinates,  produces  the  Helmholtz  equation, 


ii.  (r?l]  + 

r  dr  \  dr  J  d<j)^  dz^ 


^  ^  ^  +^  +  kln^{r,  z,  (j))j)  =  -4:-kPo5{x  -  xs)  , 


(IV.  10) 


where  ^  is  the  reference  wavenumber,  n(r,  z,  4>)  =  is  the  acoustic  index  of 

refraction,  Cq  is  the  reference  sound  speed,  and  c(r,  z,  <j))  is  the  acoustic  sound  speed. 
Most  of  the  properties  of  the  environment  are  provided  by  c(r,  z,  4>).  Regarding  the 
source  function,  it  is  assumed  to  be  a  point  source  located  at  a  distance  of  r  =  Om 
and  depth  z  =  zs  and  with  reference  source  level  Po.  This  reference  soiurce  level  is 
the  pressure  ampHtude  at  a  reference  distance  of  Ro  =  Im,  and 


<5((x))  = 


1 

27rr 


5{z-  zs)S{r) 


(rv.ii) 


is  the  Dirac-delta  function  defining  the  point  source  contribution. 
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For  simplification  of  the  Helmholtz  equation,  we  assinne  that  the  cyhndrical 
spreading  dominates  the  propagation,  and  the  magnitude  of  the  pressure  vector  will 
be  given  by 

Pir,z)  =  ^u(r,z).  (IV.  12) 

Substituting  Eq.  (IV.  12)  into  the  Helmholtz  equation  and  neglecting  the  source  term, 
we  end  up  with  the  following  equation 


d'^u  1  d%  d^u  ,2(2,  1 


u  =  0. 


(IV.  13) 


In  Eq.  (IV.  13),  we  observe  that  both  the  final  term  and  the  second  term  ,  which  intro¬ 
duces  azimuthal  couphng  between  different  radials,  are  proportional  to  — .  Therefore 

invoking  the  uncoupled  azimuth  approximation,  these  terms  will  be  neglected  in  this 
analysis. 

Defining  the  operator  notation 


P  —  — 
~  dr 


(IV.  14) 


where 


Qop  —  (//  -l-  £  +  1)^*^^  , 


2  1  1 


(IV.  15) 


(IV.  16) 


the  homogeneous  form  of  Eq.  (IV.  10)  then  becomes 

{Plp  +  KQlp)u  =  0  .  (IV.  17) 

Prop6r  factorization  of  the  outward  propagation  field  is  obtained  by  defining 

u  =  (IV.  18) 

Because  of  the  small  range  dependence  of  the  environment  it  is  also  assumed  that  the 
commutator  [Pop,Qcyp]  is  neghgible.  The  outgoing  wave  then  satisfies  [Ref.  14,  15]. 


(IV.19) 
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The  physical  meaning  of  the  Eq.  (IV.  19)  is  the  representation  acoustic  energy 
in  the  waveguide,  particularly  the  forward  propagating  acoustic  energy  when  the 
backscattered  energy  is  not  significant.  In  other  words  Eq.  (IV.19)  is  considered 
to  be  the  basic  equation  upon  which  all  one-way  PE  models  are  “built”.  The  next 
step  is  to  create  a  method  for  generating  solutions  to  this  equation  by  developing 
approximations  to  the  pseudo-differential  operator  Qop. 

For  the  creation  of  the  numerical  algorithm  which  will  solve  the  PE,  we  observe 
that  the  acoustic  field  can  be  divided  into  a  slowly  modulating  envelope  function 
and  a  phase  term  which  oscillates  at  the  acoustic  frequency.  The  PE  field  function, 
'ip{r,  z,  <f>),  is  defined  as 

^  =  V’e*''””  (IV.20) 

or,  in  terms  of  the  acoustic  pressure, 

p(r,z,4>)  =  (IV.21) 

V  r 

The  last  equation  is  scaled  in  such  a  way  that  at  r  =  Ro,  \ip\  =  1  and  |p|  =  Pq.  Com¬ 
bining  Eq.  (IV.21)  with  the  Helmholtz  equation  will  provide  us  with  wave  equation 
for  the  PE  field  function, 

^  =  —iko'tlJ  +  ikoQopfp  =  -ikoHop'ip,  (IV.22) 

or 

where 

=  (IV.23) 

is  a  Hamiltonian-like  operator  which  defines  the  evolution  of  the  PE  field  function  in 
range. 

In  Eq.  (IV.22),  the  function  is  a  vector  (in  z)  in  Hilbert  space.  The  rela¬ 
tionship  between  the  values  of  ip  at  different  ranges  can  now  be  expressed  as 

ip{r  +  Ar)  =  ^{r)ip{r).  (IV.24) 
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To  propagate  the  solution  out  in  range  requires  a  representation  of  the  propagator 
$(r).  The  MMPE  model  uses  the  spht-step  Fourier  (PE/SSF)  method  to  compute 
PE  solutions. 

In  the  SSF  algorithm  the  operators  H^p  and  Qop  are  a  combination  of  scalar 
and  differential  operators.  For  the  appropriate  functionality  of  the  SSF  algorithm 
the  different  terms  within  H^p  should  be  separated,  requring  an  approximation  to  the 
square-root  operator.  In  the  MMPE  model,  the  wide-angle  PE  (WAPE)  approxima¬ 
tion  of  Thompson  and  Chapman  is  employed  [Ref.  14,  15],  such  that 


Hop  ^  Top  +  Uop 


(IV.25) 


where 


and 


Top  =  l- 


1  + 


k^dz2 


1/2 


(IV.26) 


Uop  =  (IV.27) 

In  this  form,  the  differential  operator  has  been  separated  from  the  index  of  refraction 
term  as  required  for  implementation  with  the  SSF  technique.  In  z-space,  the  operator 
Uop  is  a.  simple  scalar  multiplication  operator,  in  other  words  a  diagonal  matrix, 
but  the  operator  Top  is  not  a  diagonal  matrix,  so  different  depth  eigenfunctions  are 
coupled.  However,  in  vertical  wavenumber  space,  the  corresponding  operator  Top  is 
diagonal. 

The  MMPE  propagator  function  is  then  based  on  the  “centered  step”  scheme. 


_  g-iko^Uop{r+Ar)^-ik„ArT,p^-iko^Uop(r) 


(IV.28) 


where  error  analysis  shows  that  this  scheme  provides  third  order  accuracy  in  Ar,  and 
is  the  method  used  in  the  MMPE  implementation.  The  general  algorithm  behind  the 
PE/SSF  implementation  is  then  as  follows.  The  PE  field  function  ^  is  specified  at 
some  range  r  in  the  z-domain.  A  multiplication  of  the  z-space  operator 
defined  at  the  beginning  of  the  range-step  is  apphed.  A  transformation  is  then  made 
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to  the  fc^-domain  followed  by  a  multiplication  of  the  fc^-space  operator  .  The 

result  is  then  transformed  again  to  the  ^-domain  followed  by  a  multiplication  of  the 
z-space  operator  defined  at  the  end  of  the  range-step.  The  final  result 

is  the  field  function  at  r  +  Ar.  The  discrete  fast  Fourier  transform  (FFT)  subroutine 
employed  in  the  numerical  code  assumes  the  convention 

i;{z)  =  FFTi-iPik,))  (IV.29) 


and 

^(k,)  =  IFFTi^iz)).  (IV.30) 


Therefore,  the  PE/SSF  implementation  can  be  represented  by 

■>p{r+Ar,z)  =  x  IFFT  ^  ^(r-,z)]}  , 

(IV.31) 


where,  in  fc^-space. 


(IV.32) 


2.  Data  Modulation-Demodulation 

The  MMPE  model  is  a  broadband,  full- wave  acoustic  propagation  model  based 
on  the  parabohc  approximation  to  the  Helmholtz  equation.  After  defining  all  the  nec¬ 
essary  parameters  for  the  entire  environment  (such  as  sound  speed  profile,  bathymetry 
of  the  water /bottom  interface,  the  acoustic  parameters  of  the  bottom,  the  soinrce  iu- 
formation),  this  model  computes  the  complex  pressure,  transmission  loss,  and  arrival 
time  structure  for  any  point  on  a  computational  grid  in  range  and  depth. 

Broadband  results  are  obtained  by  running  the  MMPE  model  for  all  discrete 
frequencies  in  the  chosen  bandwidth.  The  travel  time  results  are  then  realized  by  per¬ 
forming  a  Fourier  synthesis  and  a  necessary  multiphcation  by  some  window  function, 
S{f).  The  complex  arrival  structure  of  the  pressure  field  can  then  be  written 

/oo  P 

S(/Wr,  2,  ^  A  S{f)e-^^9(r,  f)e-^'‘df. 

(IV.33) 
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By  defining  the  reduced  time  T  =  t  -  the  phase  factor  can  be 

absorbed,  such  that 

P  P  z*®® 

P{r,  z,  T)  =  S{f)^{r,  z,  (IV.34) 

Note  that  the  use  of  reduced  time  does  not  influence  the  autocorrelation  function. 

For  all  the  environments  tested  the  same  window  function  was  used.  The 
modeled  transmitted  pulse  had  a  center  frequency  of  400  Hz  with  a  bandwidth  of 
128  Hz  and  with  256  frequency  bins  computed  in  the  DFT.  This  creates  a  signal 
of  duration  2  sec  in  the  time-domain.  In  order  to  obtain  the  real  pressure  for  each 
sonobuoy  (receiver)  while  maintaining  the  initial  assumption  that  all  the  receivers 
were  receiving  the  signals  of  the  same  duration,  we  developed  the  following  procedure. 

•  First  the  frequency-domain  complex  pressure  vector,  corresponding  to  each 
receiver,  was  multiplied  by  a  “Hanwin”  factor.  This  “Hanwin”  factor  is  a 
Harming  window  with  length  equal  of  the  number  of  the  frequency  bins. 

•  In  order  to  add  the  carrier  in  the  baseband  signal,  we  had  to  append  to  both 
sides  of  the  frequency-domain  vector  a  number  of  zeros,  creating  a  total  length 
of  1024  bins.  This  will  assist  us  with  the  frequency  resolution  that  we  are  going 
to  work  with  in  order  not  to  lose  any  information  from  the  original  transient. 
In  this  way,  that  part  of  spectrum  corresponding  to  positive  frequencies  is 
formed. 

•  The  part  of  the  spectrum  corresponding  to  negative  frequencies  was  formed 
by  taking  the  conjugate  reverse  of  the  original  vector  and  appending  it  to  the 
original  vector.  This  produced  a  complete  firequency-domain  representation  of 
the  measured  signal  over  2048  firequency  bins.  The  corresponding  sampling 
frequency  on  the  receiver  was  then  512  Hz. 

•  Performing  an  FFT  on  this  complex  frequency-domain  signal  produced  a  real 
time-domain  signal  that  represented  the  measurement  on  the  receiver.  Due  to 
a  small  imaginary  component  due  to  the  numerical  processing,  the  real  part 
of  this  signal  was  extracted  for  use  in  the  TDOA  algorithms. 

•  Before  adding  noise  to  the  transients,  all  transients  must  be  in  absolute  time 
and  not  in  reduced  time,  as  in  the  previous  step.  For  this  procedure,  a  number 
of  zeros  are  pre-appended,  which  for  each  pair  of  receivers  is  equal  to  the 
fraction  of  the  difference  of  the  two  distances  (of  the  2  receivers)  divided  by 
the  factor  dt  =  1/1024. 


36 


3.  MMPE  Environments 

Before  proceeding  with  the  simulation  results,  let  us  describe  the  four  test 
cases  that  were  used  with  the  MMPE  model.  This  will  provide  a  common  reference 
not  only  for  the  TDOA  simulation  results  but  also  for  the  locahzation  and  tracking 
problems  presented  in  the  next  chapters.  The  four  environments  are: 

1.  Flat  Bottom 

2.  Sound.  Channel 

3.  Shelf  Break 

4.  Internal  Waves 

The  first  two  of  these  are  range-independent  (RI)  while  the  last  two  are  range- 
dependent  (RD).  The  main  difference  between  these  two  categories  is  that  in  the 
range-independent  case  the  received  signal  at  each  sonobuoy  depends  only  on  the 
distance  between  the  source  and  the  receiver.  In  other  words,  the  sound  speed  profile 
and  bottom  are  the  same  at  every  point  in  the  area  of  interest.  In  the  range-dependent 
case  the  sound  speed  profile  and/or  bottom  changes  according  to  range  and  direction 
from  the  source.  This  makes  it  more  difficult  to  simulate  since  one  has  to  be  aware  in 
more  detail  about  the  environmental  conditions.  Range-dependent  simulations  pro¬ 
vide  more  realistic  data  consistent  with  the  complicated  conditions  that  occur  in  the 
ocean. 

а.  Flat  Bottom 

This  environment  is  defined  by  a  flat  ocean  bottom  and  isospeed  water 
colrunn.  In  our  case  we  have  regularly  spaced  values  of  bottom  parameters.  At  each 
range  location,  the  bottom  is  assumed  to  be  a  homogeneous  half-space  below  that 
position.  The  specific  values  of  the  properties  for  the  case  are  listed  in  Table  II: 

б.  Sound  Channel 

This  range-independent  case  is  similar  to  the  previous  one  with  the 
following  major  differences.  First,  the  bottom  bathymetry  is  varying  and  the  sound 
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Bathymetry 

flat  (300  m  depth) 

Soimd  Speed  Profiles 

isospeed  (1500  m/s) 

Water  Density 

1  g/cm^ 

Water  Attenuation 

0.0  dB/A 

Bottom  Sound-Speed 

1700  m/s 

Bottom  Density 

1.6  g/cm^ 

Compressional  Bottom 
Attenuation 

0.1  dB/km/Hz 

Shear 

None 

Table  II.  Environmental  Properties  for  RI  Test  case:  “Flat  Bottom”. 


speed  profile,  which  is  the  same  for  the  entire  computational  grid,  is  defined  by  the 
following  equations: 


c=  < 


1515  +  0.016z 


1490 


1  +  0.25  sw*' )  + 


when  Z  G  [Zg^urf^  ^duct] 

,  (IV.35) 

when  z  G  [z^^ict,  z^ax]  ? 


where  Zg^rf  —  Om  is  the  surface  depth,  Zdua  —  75m  is  the  depth  of  the  surface  duct, 
Zmax  =  500m  is  the  maximum  depth  and  Zaxis  =  200m  is  the  depth  of  the  axis  of  the 
sound  channel.  The  soxmd  speed  profile  is  shown  in  Fig.  7.  The  specific  values  of  the 
properties  for  the  case  are  listed  below  in  Table  III. 


Bathymetry 

varying 

Sound  Speed  Profiles 

varying 

Water  Density 

i  g/cm^ 

Water  Attenuation 

0.0  dB/A 

Bottom  Sound-Speed 

1700  m/s 

Bottom  Density 

1.6  g/cm^ 

Compressional  Bottom 
Attenuation 

0.1  dB/km/Hz 

Shear 

None 

Table  III.  Environmental  Properties  for  RI  Test  case:  “Sovmd  Charmel”. 
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Figure  7.  Sound  Speed  Profile  for  RI  Test  case:  Sound  Channel, 
c.  Shelf  Break 

The  next  case  has  to  do  with  a  range-dependent  environment  similar 
to  that  foimd  in  the  regions  of  continental  shelf  breaks.  This  environment  is  defined 
by  both  a  varying  water  column  sound  speed  and  bottom  bathymetry  with  homoge¬ 
neous  bottom  acoustical  properties.  This  environment  is  loosely  based  on  shelf  break 
front  structures  observed  during  the  ONR  Primer  experiment  [Ref.  31].  The  general 
properties  of  this  test  case  are  defined  in  Table  IV. 


Bathymetry 

varying 

Sound  Speed  Profiles 

varying 

Water  Density 

1  g/cm^ 

Water  Attenuation 

0.0  dB/A 

Bottom  Sound-Speed 

1700  m/s 

Bottom  Density 

1.5  g/cm^ 

Compressional  Bottom 
Attenuation 

0.1  dB/km/Hz 

Shear 

None 

Table  IV.  Environmental  Properties  for  RD  Test  case:  “Shelf  Break” . 


39 


The  vertical  sound  speed  profile  and  bathymetry  for  this  case  are  dis¬ 
played  in  Fig.  8(a).  This  characteristic  is  extended  infinitely  in  the  direction  perpen¬ 
dicular  to  the  page  providing  a  3-D  shelf  break  environment.  The  locations  of  the 
source  and  two  of  the  receivers  are  shown  in  Fig.  8(b),  which  is  the  same  environment 
as  Fig.  8(a)  but  viewed  from  the  top.  Isobaths  are  also  drawn  to  provide  a  better 
perspective  on  the  relative  positions  of  the  source  and  receivers. 


0  5  10  16  20 

Range  (km) 


(a) 

Isobaths  (m) 


400  350  300  250  200  150  100  100 


Range  (km) 

(b) 

Figure  8.  Sound  speed  profile:  (a)  Shelf  break  environment  vertical  view,  (b)  Shelf 
break  environment  horizontal  view  with  sample  locations  of  source  and  receivers. 
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d.  Internal  Waves 

This  environment  is  defined  by  varying  water  colmnn  sound  speed  fea¬ 
tures  and  a  fiat,  homogeneous  bottom.  More  specifically  this  test  case  is  the  combi¬ 
nation  of  a  simple,  sinusoidally  varying  perturbation  with  a  soliton  wavefield,  and  is 
loosely  based  on  the  work  of  Tielburger,  Finette,  and  Wolf  [Ref.  32].  The  aforemen¬ 
tioned  environment  is  characterized  by  the  following  properties  as  shown  in  Table  V. 


Bathymetry 

Flat  bottom  (zb  =  200  m  depth) 

Sound  Speed  Profiles 

varying 

Water  Density 

1  g/cm^ 

Water  Attenuation 

0.0  dB/A 

Bottom  Sound-Speed 

1700  m/s 

Bottom  Density 

1.5  g/cm^ 

Compressional  Bottom 
Attenuation 

0.1  dB/km/Hz 

Shear 

None 

Table  V.  Environmental  Properties  for  RD  Test  case;  “Internal  Waves” . 


The  background  (range-independent)  sound  speed  profile  used  is  the 
same  as  defined  in  Eq.  (IV.35).  The  perturbations  to  this  background  are  a  com¬ 
bination  of  a  sum  of  simple  sinusoids  and  a  train  of  soliton  waves.  The  sinusoidal 
perturbations  are  defined  by 

dcsin{z,r)  =  C  ^  {cos  (ii'(z)r)}  (IV.36) 

i=i 


where  K{i)  =  [(2ooom)-(i^i)(3bbm)i>  ^  ~  25m,  and  C  is  defined  such  that  the  maximvun 
value  for  dcsin  is  7.bm/s. 

The  sohton  perturbation  is  defined  by 


dcsol{z,r)  =  C  X!  I 

i=i  I 


seek 


V  m  . 


(IV.37) 
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where  B  =  25m,  A(i)  =  10e(-°-3(*-i)),  D{t)  =  Rii)  =  Rii  -  1)  -  (7  -  ^500 

for  z  =  2,  •  •  • ,  6,  i2(l)  —  14000m,  and  C  is  defined  such  that  the  maximum  value  for 
dc  is  12.5m/s. 

The  combination  of  the  above  perturbations  is  then  simply  defined  by 

dc{z,  r)  =  dcsol{z,  r)  +  dcsin{z,  r)  (IV.38) 

and  the  result  for  the  total  sound  speed  structure  is  depicted  in  Fig.  9.  Note  that 
the  solitons  are  modeled  as  propagating  only  in  the  plane  of  the  page  while  the 
internal  wave  sinusoids  are  considered  randomly  oriented  in  all  directions.  For  the 
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Figure  9.  Sound  speed  profile:  Combination  of  Sinusoidal  and  Soliton  Perturbations. 

range-independent  environments  the  source  depth  was  between  140m  and  160m  and 
the  locations  of  the  receivers  were  at  ranges  between  2km  and  8km  and  at  depths 
between  40m  and  60m.  For  the  Shelf  Break  environment  the  source  depth  was  at 
250m  and  the  receivers  were  at  a  depth  of  100m  at  ranges  between  2km  and  5km, 
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while  for  the  Internal  Waves  environment  the  source  depth  was  at  150m  and  the 
receivers  were  at  a  depth  of  50m  and  at  the  same  ranges  as  in  the  Shelf  Break.  Better 
insight  to  the  scenarios  is  given  in  the  next  chapter,  where  the  simulation  results  for 
the  TDOA  are  presented. 
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V.  TDOA  SIMULATION  RESULTS 


A.  SYNTHETIC  DATA 

In  this  section  we  present  simulation  results  using  synthetic  data  in  the  TDOA 
estimation  problem.  Our  purpose  is  to  distinguish  how  the  two  famihes  of  methods 
(classical  -  subspace)  perform,  what  their  basic  differences  are,  and  finally  which 
method(s)  produce  the  best  overall  results.  Since  the  number  of  combinations  using 
different  signals,  kind  and  amount  of  noise,  and  different  characteristics  of  the  signals 
and  of  the  methods  is  huge,  we  have  tried  to  provide  some  typical  cases  in  order  to 
provide  insight  into  the  behavior  of  the  methods.  We  will  try  to  explain  how  the 
different  factors  attribute  to  the  final  result,  which  is  the  estimation  of  TDOA. 

In  the  first  subsection  of  the  synthetic  data  we  present  cases  from  “Expo¬ 
nential”,  “Sinusoidal”,  “Damped  Sinusoidal”,  “Chirp”  and  finally  “Damped  Chirp” 
signals  at  different  noise  levels  (SNR  15,  10,  5  dB).  Our  intent  is  to  show  the  be¬ 
havior  of  the  methods  when  the  transient  duration  changes  (short  to  long)  and  how 
successful  they  are  when  the  desired  TDOA  to  be  predicted  is  large  or  small.  In 
all  cases  there  are  different  realizations  of  white  noise  with  different  variance  added 
to  each  transient.  For  the  subspace  methods,  the  size  of  the  covariance  matrix  (see 
Chapter  3)  taken  was  10  and  20  samples,  and  for  the  classical  methods  the  munber  of 
segments  (see  Chapter  2)  was  either  1  or  2.  In  addition  to  the  figures,  corresponding 
tables  summarize  the  results  in  order  to  provide  a  quick  appreciation  of  each  method’s 
behavior. 

In  the  second  subsection  of  the  synthetic  data,  we  merge  both  kinds  of  methods 
(classical-subspace)  to  examine  how  the  category  of  methods  influences  the  other  and, 
of  course,  if  there  is  any  substantial  improvement  in  the  results. 

1.  Implementation  Using  Individual  Methods 

The  cases,  whose  results  are  presented  in  this  subsection  and  in  Appendix  A, 
using  synthetic  data  are  the  following: 
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•  exponential  transient  with  length  of  100s  and  actual  TDOA  50s; 

•  exponential  transient  with  length  of  20s  and  actual  TDOA  100s; 

•  exponential  transient  with  length  of  50s  and  actual  TDOA  -10s; 

•  sinusoidal  transient  with  length  of  100s  and  actual  TDOA  50s; 

•  damped  sinusoidal  transient  with  length  of  100s  and  actual  TDOA  50s; 

•  chirp  transient  with  length  of  100s  and  actual  TDOA  50s;  and 

•  damped  chirp  transient  with  length  of  100s  and  actual  TDOA  50s  . 

Each  case  was  tested  with  different  realizations  of  additive  noise  and  under  three 
different  levels  of  SNR.  The  figures  and  the  comparative  table  for  the  first  case  is 
presented  in  this  subsection.  The  remaining  case  results  are  submitted  in  Appendix 
A.  Even  though  we  have  tested  each  kind  of  transient  thoroughly  changing  various 
parameters,  apart  from  the  ones  already  defined,  we  have  observed  that  the  major 
influence  to  the  desired  result  of  TDOA  is  provoked  by  the  noise  (SNR),  length 
of  transient,  amount  of  predicted  TDOA,  size  of  covariance  matrix  (for  subspace 
methods),  number  of  segments  (for  classical  methods),  and  kind  of  transient.  Having 
all  the  above  in  mind,  we  picked  three  typical  cases  from  one  kind  of  transient  in  order 
to  see  how  the  methods  perform  under  different  values  of  the  aforementioned  factors, 
and  the  same  case  from  the  rest  of  the  transients  for  checking  possible  difficulties  in 
estimating  the  TDOA  for  a  particular  type  of  signal. 

Figures  10  through  15  show  TDOA  results  for  the  ejqponential  transient  for 
all  methods  using  SNR  values  of  15dB,  lOdB,  and  5dB.  The  estimated  values  of  time 
delay  are  fisted  in  Table  VI. 
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Figure  10.  Exponential  transient:  length  L==100s;  Actual  TDOA=50s;  SNR=15dB; 
White-Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods, 
number  of  segments=l. 
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Figure  11.  Exponential  transient:  length  L=100s;  Actual  TDOA=50s;  SNR=15dB 
White-Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods 
number  of  segments=2. 


Figiire  12.  Exponential  transient:  length  L=100s;  Actual  TDOA=50s;  SNR=10dB; 
White-Noise  (a^).  (a)  Subspace  methods,  Covariance  Size=10  (b)  Classical  methods, 
number  of  segments=l 
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Figure  13.  Exponential  transient:  length  L=100s;  Actual  TDOA=50s;  SNR=10dB; 
White-Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods, 
number  of  segments=2. 


Figure  14.  Exponential  transient:  length  L==100s;  Actual  TDOA=50s;  SNR=05dB; 
White-Noise  (cTg).  (a)  Subspace  methods,  Covariance  Size=10  (b)  Classical  methods, 
number  of  segments=l 
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Figure  15.  Exponential  transient:  length  L=100s;  Actual  TDOA=50s;  SNR=05dB; 
White-Noise  {(rl).  (a)  Subspace  methods,  Covariance  Size=20  (b)  Classical  methods, 
number  of  segments=2 


Method 

Parameters 

SNR  15  dB 

SNR  10  dB 

SNR  5  dB 

MUSIC 

Cov-Mat:  10 

50 

51 

51 

Cov-Mat:  20 

50 

49 

49 

Modified  MUSIC 

Cov-Mat:  10 

50 

51 

52 

Cov-Mat:  20 

50 

49 

50 

MIN-NORM 

Cov-Mat:  10 

50 

51 

50 

Cov-Mat:  20 

50 

49 

49 

PCLP 

Cov-Mat:  10 

50 

51 

52 

Cov-Mat:  20 

50 

49 

49 

ESPRIT 

Cov-Mat:  10 

49 

51 

50 

Cov-Mat:  20 

50 

50 

50 

Root-MUSIC 

Cov-Mat:  10 

50 

50 

49 

Cov-Mat:  20 

50 

49 

49 

SCOT 

#-Segs:  1 

50 

50 

92 

#-Segs:  2 

50 

52 

47 

PHAT 

#-Segs:  1 

50 

50 

92 

#-Segs:  2 

50 

52 

47 

Roth 

#-Segs:  1 

-299 

-502 

248 

#-Segs:  2 

50 

52 

-204 

X-Corr 

#-Segs:  1 

50 

50 

53 

#-Segs:  2 

50 

50 

47 

Table  VI.  Exponential  transient:  length  L=100s;  white  noise  =  2);  actual 
TDOA=50s. 

After  examining  the  corresponding  Figs.  10  - 15  and  Table  VI  for  the  first  case 
of  the  synthetic  data  we  can  say  that  almost  all  the  methods  performed  satisfactorily. 
The  figures  for  the  subspace  methods  in  all  cases,  Figs.  10(a)  -  15(a),  were  clear 
without  any  subsidiary  peaks,  compared  with  the  respective  figures  of  the  classical 
methods.  Figs.  10(b)  -  15(b).  For  medium  noise  level  (SNR=10dB),  Figs.  12  -  13, 
and  even  more  for  high  noise  level  (SNR=5dB),  Figs.  14  -  15,  the  subspace  methods 
gave  more  accurate  results.  The  increase  of  the  size  of  the  covariance  matrix  had  a 
slight  improvement  in  the  performance  for  the  subspace  methods.  On  the  contrary 
the  increase  of  the  number  of  segments  for  the  classical  methods  had  more  positive 
influence  on  their  performance,  especially  when  the  SNR  became  very  low  (5dB).  Of 
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the  classical  methods,  Roth  had  the  worst  performance  and  cross-correlation  had  the 
best  performance,  on  average.  Subspace  methods  with  both  approaches  had  similar 
behavior  in  all  conditions. 

In  summary,  after  reviewing  all  cases  of  the  synthetic  data  (see  Figs.  51  - 
86  and  Tables  XX  -  XXIV  in  Appendix  A)  for  the  performance  of  all  methods,  the 
following  observations  can  be  made. 

•  The  subspace  methods  uniformly  had  very  efficient  performance  compared 
with  the  classical  methods 

•  The  subspace  methods  were  more  consistent  in  the  quality  of  the  results  as 
a  group.  The  classical  methods,  on  the  other  hand,  in  some  cases  had  great 
diversity  among  them  in  the  quality  of  their  performance.  Among  the  clas¬ 
sical  methods,  Roth  had  the  poorest  performance;  and  in  general  the  cross 
correlation  was  more  consistent  in  accmacy  than  the  other  methods. 

•  The  plots  of  the  subspace  methods  were  more  distinct  and  it  was  easier  to  iden¬ 
tify  the  peak  corresponding  to  the  correct  Time  Difference  of  Arrival  (TDOA). 

•  In  cases  where  both  methods  provided  the  expected  results  the  classical  meth¬ 
ods  had  more  acciurate  peaks  than  the  subspace  methods.  However  both  meth¬ 
ods  provided  results  in  which  accuracy  was  completely  acceptable. 

•  In  order  to  perform  with  multiple  segments  the  classical  methods  needed  to 
have  appropriate  length  of  input  data  (signal  plus  noise).  Otherwise  phe¬ 
nomena  of  aliasing  appeared  especially  in  cases  with  a  large  “value”  of  TDOA 
between  the  two  data  sets.  This  drawback  did  not  affect  the  subspace  methods 
at  all  since  they  do  not  deal  with  segmentation. 

•  Methods  like  root  MUSIC  or  ESPRIT  can  produce  direct  numerical  estimates 
for  the  time  delay  without  the  need  to  search  for  a  peak. 

•  Subspace  methods  were  more  “resistant”  to  noise  compared  to  the  classical 
methods. 

•  For  short  duration  of  transients  the  performance  of  the  methods  decreased  in 
accuracy,  especially  when  it  was  combined  with  low  SNR.  In  this  case  classical 
methods  gave  totally  wrong  results  compared  with  the  subspace  methods, 
whose  performance  was  inside  the  permissible  hmits  of  error  (Table  XX). 

•  For  the  case  of  a  small  value  of  predicted  TDOA  (Table  XXI)  the  subspace 
methods  performed  more  than  satisfactorily  even  under  high  “noise  condi¬ 
tions”.  The  performance  of  the  classical  methods  was  similar  with  that  of 
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the  subspace  methods  for  high  SNR,  but  they  became  quite  inferior  when  the 
noise  was  increased. 

•  In  general  the  increase  in  the  size  of  the  covariance  matrix  in  subspace  methods 
improved  their  performance,  especially  when  the  SNR  was  small.  Better  results 
were  achieved  when  the  size  of  the  covariance  matrix  was  greater  than  20,  but 
increased  the  total  niunber  of  computations  something  that  made  the  subspace 
methods  to  delay  enough  in  order  to  provide  us  with  the  desired  result. 

For  the  cases  where  all  the  remaining  factors  are  the  same  (except  for  the  type  of 
transient).  Tables  XXII,  XXIII,  XXIV,  and  Figs.  10  -  15  and  63  -  86,  we  can  say 
that  all  the  methods  have  almost  the  same  behavior  with  different  transients  as  with 
the  exponential  transient.  Only  the  sinusoidal  transient  had  noticeably  worse  perfor¬ 
mance  than  the  others  and  was  influenced  more  with  high  noise  level  than  the  rest  of 
the  transients. 

2.  Implementation  Using  Combination  of  Methods 

Before  we  continue  with  the  MMPE  data,  it  is  worth  examining  if  the  “com¬ 
bination”  of  the  two  families  of  methods,  subspace  and  classical,  perform  better  than 
each  one  separately.  In  our  case  the  word  “combination”  means  that  before  we  apply 
the  subspace  methods  in  the  frequency  domain  (cross-spectrum),  the  desired  weight¬ 
ing  factor  will  be  implemented  according  to  the  processor  we  would  like  to  use.  The 
main  idea  was  to  examine  if  the  subspace  methods,  with  the  assistance  of  the  weight¬ 
ing  factor,  could  provide  us  with  more  accurate  results  (peaks).  Several  test  cases 
were  conducted  and  it  was  determined  that  the  peaks  were  more  accurate  but  the 
results  were  not  better  than  without  implementing  the  combination.  As  a  general 
remark,  each  weighting  factor  influenced  the  performance  of  the  subspace  methods 
differently  in  such  a  way  that  the  subspace  methods  lost  their  good  consistency  in 
results.  Figures  16-17  and  the  corresponding  Table  VII  give  a  sample  of  the  “com¬ 
bination”  performance. 
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Method 

SCOT 

PHAT 

ROTH 

CROSS 

CORRELATION 

NONE 

MUSIC 

95 

95 

-879 

103 

103 

Modified-MUSIC 

96 

96 

-600 

103 

103 

MIN-NORM 

90 

90 

-620 

103 

103 

PCLP 

91 

91 

-871 

103 

103 

ESPRIT 

95 

95 

-728 

103 

103 

Root-MUSIC 

96 

96 

-864 

103 

103 

Table  VII.  Damped  Sinusoidal  transient:  length  L=400s;  white  noise  (cr^  =  2); 
SNR=10dB;  actual  TDOA=100s. 
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Figure  16.  Damped  Sinusoidal  transient:  Actual  TDOA=100s;  SNR=10dB;  White- 
Noise  ((Tq)  Covariance  Size=10,  number  of  segments=l.  (a)  Subspace  methods  with 
SCOT,  (b)  Subspace  methods  with  PHAT. 
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Figure  17.  Damped  Sinusoidal  transient:  Actual  TDOA=100s;  SNR=10dB;  White- 
Noise  ((Jo)  Covariance  Size=10,  number  of  segments=l.  (a)  Subspace  methods  with 
ROTH,  (b)  Subspace  methods  with  Cross  Correlation. 
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B.  MMPE  DATA 

In  this  section  results  from  data  generated  by  the  MMPE  propagation  model 
are  presented.  The  particular  environments  that  were  used  to  generate  these  data 
are  the  ones  presented  in  Chapter  4.  These  simulations  were  conducted  to  verify 
the  observations  from  the  synthetic  data  and  to  test  the  methods  with  broadband 
transients  under  more  complicated  environmental  conditions.  For  each  environment, 
we  provide  a  Transmission  Loss  (TL)  plot  relative  to  the  depth  and  the  frequency 
bandwidth,  the  waveforms  of  the  source  and  of  the  two  receivers  in  the  time  domain, 
and  the  results  of  all  TDOA  methods.  It  will  be  seen  that  in  all  cases  the  signals 
in  each  location  have  totally  different  shape  and  magnitude  due  to  the  multipath 
structure  of  the  waveguide.  As  mentioned  before,  the  environments  are  divided  into 
two  categories:  range-independent  (the  received  signal  depends  only  on  the  distance 
between  the  source  and  the  receiver  [Flat  -  Sound  Channel  cases],  and  range- dependent 
(the  sound  speed  profile  changes  according  to  range  and  direction  from  the  source 
[Shelf  Break  -  Internal  Waves  cases].  Two  cases  are  considered  for  each  environment. 
For  the  range-independent  environments  the  receivers  are  chosen  to  be  at  different 
ranges  from  the  source.  For  the  range-dependent  environments  the  receivers  are 
selected  to  be  at  the  same  ranges  but  in  different  relative  directions  from  the  source. 
In  this  way  the  difference  between  the  two  categories  will  be  more  evident.  For  more 
detailed  description  of  the  MMPE  environments  see  Chapter  4. 

1.  Flat  Bottom 

Two  cases  are  considered  for  this  environment.  The  receivers  are  at  distances 
of  2km  and  3km,  and  2km  and  5km  for  each  case,  respectively.  In  all  cases  different 
realizations  of  white  noise  were  added  to  each  transient  with  SNR=10dB.  In  Figs.  18 
and  19  are  shown  the  Transmission  Loss  (TL)  plot  of  depth  versus  frequency  and  the 
waveforms  at  the  source  and  at  the  two  receivers  in  the  time  domain  for  one  of  the 
two  cases.  The  results  for  both  cases  are  shown  in  Figs.  20  and  21  and  summarized 
in  Table  VIII.  Prom  these  results,  it  appears  that  all  methods  perform  in  a  similar 
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manner  to  that  in  the  section  with  the  synthetic  data.  All  of  them  gave  results  very 
close  to  the  actual  value  for  both  cases  with  shghtly  better  accuracy  by  the  subspace 
methods.  The  difference  in  the  figures  between  the  two  group  of  methods  is  still 
the  same.  Subspace  methods  have  more  clear  plots  and  it  is  easier  to  pick  up  the 
correct  peak  from  their  graphs  than  the  corresponding  ones  from  classical  methods 
(subsidiary  peaks.)  In  this  environment  the  root-finding  subspace  methods  (ESPRIT 
-  rootMUSIC)  had  the  best  performance  of  all  the  methods  and  Roth  gave  the  worst 
results  compared  with  the  rest. 


Method 

Estimated  TDOA 

Estimated  TDOA 

Case  1 

Case  2 

ACTUAL 

0.6936 

2.1209 

MUSIC 

0.6953 

2.1250 

Modified-MUSIC 

0.6641 

2.0879 

MIN-NORM 

0.6943 

2.1240 

PCLP 

0.6943 

2.1240 

ESPRIT 

0.6934 

2.1221 

Root-MUSIC 

0.6943 

2.1230 

SCOT 

0.7324 

2.1182 

PHAT 

0.7324 

2.1182 

Roth 

0.5713 

2.1328 

X-Corr 

0.7285 

2.1133 

Table  VIII.  TDOA  results  for  Flat  Bottom  isospeed  case. 
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Figure  18.  Transmission  Loss  for  FLAT  Bottom  environment. 
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Figure  19.  Time  domain  waveforms  for  FLAT  Bottom  environment,  (a)  Source,  (b) 
Receiver  1.  (c)  Receiver  2. 
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T<sec)  T(seo) 

(b) 

Figure  21.  FLAT  Bottom  environment  Subcase  2:  Actual  TDOA=2.1209s; 

White-Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=10.  (b)  Clas¬ 
sical  methods,  number  of  segments=l. 


2.  Sound  Channel 

Two  cases  are  considered  for  this  environment.  In  the  first  case  the  receivers 
are  at  distances  2km  and  4km,  while  in  the  second  case  they  are  at  2km  and  8km.  In 
all  cases  a  different  realization  of  white  noise  was  added  to  each  transient  to  obtain 
a  SNR  of  lOdB.  Figures  22  and  23  show  the  Transmission  Loss  (TL)  plot  of  depth 
versus  frequency  and  the  waveforms  of  the  source  and  of  the  two  receivers  in  the 
time  domain  for  one  of  the  two  cases.  The  results  for  both  cases  are  shown  in  Figs. 
24  and  25  and  summarized  in  Table  IX.  The  results  show  us  a  similar  picture  for 
the  performance  of  the  methods.  The  only  difference  is  that  in  this  environment  the 
“gap”  between  the  estimated  values  of  TDOA  from  the  classical  methods  and  the 
actual  predicted  value  becomes  larger. 


Method 

Estimated  TDOA 

Estimated  TDOA 

Case  1 

Ceise  2 

ACTUAL 

1.5830 

4-5929 

MUSIC 

1.5830 

4.5967 

Modified-MUSIC 

1.5781 

4.5967 

MIN-NORM 

1.5830 

4.5967 

PCLP 

1.5830 

4.5967 

ESPRIT 

1.5820 

4.5967 

Root-MUSIC 

1.5820 

4.5967 

SCOT 

1.5771 

4.5527 

PHAT 

1.5571 

4.5527 

Rolh 

1.5596 

4.5527 

X-Corr 

1.5732 

4.5527 

Table  IX.  TDOA  results  for  SOUND  CHANNEL  case. 
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Figure  22.  Transmission  Loss  for  SOUND  CHANNEL  environment. 
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Figure  23.  Time  domaiu  waveforms  for  SOUND  CHANNEL  environment,  (a)  Source, 
(b)  Receiver  1.  (c)  Receiver  2. 
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Figure  24.  SOUND  CHANNEL  environment  Subcase  1:  Actual  TDOA=1.5830s; 
SNR=10dB;  White-Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=10.  (b) 
Classical  methods,  number  of  segments=l. 
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Figure  25.  SOUND  CHANNEL  environment  Subcase  2:  Actual  TDOA=4.5967s; 
SNR=10dB;  White-Noise  (a^).  (a)  Subspace  methods,  Covariance  Size=10.  (b) 
Classical  methods,  number  of  segments=l. 


3.  Shelf  Break 

For  this  environment  the  receivers  are  at  distances  2km  and  3km  but  in  posi¬ 
tions  relative  to  the  source  as  shown  in  Fig.  26  for  each  of  the  two  cases  considered. 
In  particular  Fig.  26  provides  us  with  a  top  view  of  the  geometry  for  the  two  cases, 
indicating  the  positions  of  the  source  and  of  the  receivers  relative  to  the  source,  noting 
also  the  corresponding  depths  and  the  respective  ranges  for  each  case.  Once  again, 
different  realizations  of  white  noise  were  added  to  each  transient  to  obtain  a  SNR 
of  lOdB.  Figures  27  and  28  show  the  Transmission  Loss  (TL)  plot  of  depth  versus 
frequency  and  the  waveforms  of  the  source  and  of  the  two  receivers  in  the  time  do¬ 
main  for  one  of  the  two  cases.  The  results  for  both  cases  are  shown  in  Figs.  29  and 
30  and  summarized  in  Table  X.  As  we  now  get  into  more  complex  environments 
(range-dependent),  it  is  obvious  that  the  graphs  of  the  classical  methods  become 
more  “noisy”  than  the  ones  in  the  range-independent  environments,  and  also  that 
they  had  an  abrupt  decrease  in  their  performance  except  for  the  Cross-correlation 
method.  Subspace  methods  still  gave  us  accurate  results  and  all  of  them  were  very 
consistent  about  the  small  diversity  of  their  estimated  values  of  TDOA. 


LEGENI> _ 

Case  1  Receivers 
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Case  2  ReceiversCj) 
Source  Depth  :  250m 

Receiver  Depth:  1 00m 


Figure  26.  Top  view  of  positions  of  source  and  receivers  for  both  cases  of  SHELF 
BREAK  environment. 
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Figure  27.  Transmission  Loss  for  SHELF  BREAK  environment. 


Method 

Estimated  TDOA 

Estimated  TDOA 

Case  1 

Case  2 

ACTUAL 

0.7173 

0.7407 

MUSIC 

0.7139 

0.7402 

Modified-MUSIC 

0.7139 

0.7490 

MIN-NORM 

0.7139 

PCLP 

0.7139 

ESPRIT 

0.7139 

0.7471 

Root-MUSIC 

0.7139 

0.7471 

SCOT 

-0.2793 

PHAT  \ 

-0.2793 

Roth 

1.2998 

1.6211 

X-Corr 

0.8154 

0.7803 

Table  X.  TDOA  results  for  SHELF  BREAK  case. 
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Figiire  28.  Time  domain  waveforms  for  SHELF  BREAK  environment,  (a)  Source, 
(b)  Receiver  1.  (c)  Receiver  2. 
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sical  methods,  number  of  segments=l. 
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4.  Internal  Waves 

For  this  environment  both  receivers  are  at  a  distance  of  3km  but  at  different 
positions  relative  to  the  source  as  shown  in  Fig  31  for  each  of  the  two  cases  considered. 
Figure  31  provides  a  top  view  of  the  geometry  for  the  two  cases,  indicating  the 
positions  of  the  somce  and  of  the  receivers  relative  to  the  source,  noting  also  the 
corresponding  depths  and  the  respective  ranges  for  each  case.  As  with  the  previous 
environments,  different  realizations  of  white  noise  were  added  to  each  transient  to 
obtain  a  SNR  of  lOdB.  Figures  32  and  33  show  the  Transmission  Loss  (TL)  plot  of 
depth  versus  frequency  and  the  waveforms  at  the  source  and  at  the  two  receivers  in 
the  time  domain  for  one  of  the  two  cases.  The  results  for  both  cases  are  shown  in  Figs. 
34  and  35  and  summarized  in  Table  XI.  Our  observations  for  this  environment  are 
almost  the  same  as  with  the  Shelf  Break.  The  only  difference  is  that  even  the  subspace 
methods  give  slightly  worse  results  compared  with  the  actual  value  of  TDOA.  This 
was  expected  since  the  Internal  Waves  was  the  most  complex  environment  used  in 
this  thesis.  These  results  are  still  within  the  permissible  limits  of  error. 


Case  1  Receivers  (2^ 

Soiirce 

O 

Case  2  Receivers  <© 

Source  Depth.  : 

150m 

Receiver  Depth: 

50m 

Figure  31.  Top  view  of  positions  of  source  and  receivers  for  both  cases  of  INTERNAL 
WAVES  environment. 
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Figure  32.  Transmission  Loss  for  INTERNAL  WAVES  environment. 


Method 

Estimated  TDOA 

Estimated  TDOA 

Case  1 

Case  2 

ACTUAL 

0.7770 

0.2109 

MUSIC 

0.7734 

0.2090 

Modified-MUSIC 

0.8311 

0.1885 

MIN-NORM 

0.7822 

0.2061 

PCLP 

0.7793 

ESPRIT 

0.7959 

RooUMUSIC 

0.7930 

0.2041 

SCOT 

0.1289 

PHAT 

0.1289 

Roth 

2.5938 

X-Corr 

0.8574  1 

0.1357 

Table  XL  TDOA  results  for  INTERNAL  WAVES  case. 
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(c) 

Figure  33.  Time  domain  waveforms  for  INTERNAL  WAVES  environment,  (a) 
Source,  (b)  Receiver  1.  (c)  Receiver  2. 
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MInNorm:  0.78223 


M  Music:  0.83105 


Figure  34.  INTERNAL  WAVES  environment  Subcase  1;  Actual  TDOA=0. 7770s; 
SNR=10dB;  White-Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=10.  (b) 

Classical  methods,  number  of  segments=l. 


Figure  35.  INTERNAL  WAVES  environment  Subcase  2:  Actual  TDOA=0.2109s; 
SNR=10dB;  White-Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=10.  (b) 
Classical  methods,  number  of  segments=l. 


Taking  into  consideration  the  figures  and  the  tables  with  the  TDOA  results 
from  aU  the  environments,  we  observe  that  aU  the  methods  performed  quite  accurately. 
It  should  be  noted  that  the  subspace  methods  were  more  consistent  than  the  classical 
ones.  This  can  be  verified  from  the  range-dependent  environments,  where  all  the 
subspace  methods  gave  results  more  close  to  the  actual  value  and  they  had  only  a  small 
dispersion  among  them.  On  the  contrary,  the  classical  methods  had  slightly  inferior 
performance  in  the  range-independent  environments  and  much  worse  in  the  range- 
dependent  ones,  with  a  wider  dispersion  of  values.  This  difference  in  the  performance 
among  all  the  methods  between  the  two  categories  of  the  environments  also  appears 
in  the  results  of  the  localization  problem  in  the  next  chapter. 
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VI.  THE  LOCALIZATION  PROBLEM 


In  this  chapter,  we  use  the  TDOAs  between  the  receivers,  estimated  by  the 
various  methods,  in  order  to  estimate  the  location  of  the  target.  For  this  purpose  a 
homonymous  algorithm  will  be  implemented  whose  main  characteristic  is  to  provide 
us  the  passive  assessment  of  target  position  by  using  a  corresponding  set  of  equations, 
called  Time  of  Arrival  (TOA)  equations.  The  TOA  equations  try  to  correlate  the 
positions  of  the  receivers  with  the  position  of  the  target  using  the  time  that  a  transient 
needs  to  travel  from  one  location  to  the  other.  In  other  words,  these  equations  try 
to  find  a  correspondence  between  the  signal  arrival  times  and  the  target  range.  By 
solving  the  TOA  equations  for  only  two  receivers,  we  are  not  going  to  get  one  specific 
location  of  the  target,  but  an  infinite  number  of  possible  locations  since  this  solution 
corresponds  to  a  bearing  from  this  peirticular  pair  of  receivers.  Under  these  conditions, 
in  order  to  eliminate  this  ambiguity  we  need  to  solve  the  TOA  equations  for  all  the 
receivers  at  the  same  time,  so  the  desired  target  position  wiU  be  the  result  of  the 
crossing  of  all  the  bearings.  This  expectation  is  quite  natural  since,  as  far  as  we 
know,  only  one  solution  (target  position)  can  account  simultaneously  for  the  travel 
time  of  the  transient  to  each  receiver.  The  algorithm  that  is  implemented  is  simpler 
than  most  since  it  transforms  the  original  nonlinear  TOA  equations  into  a  set  of  linear 
equations.  The  inputs  (receiver  position  and  signal’s  arrival  time  for  each  receiver) 
correspond  to  a  specific  set  of  outputs  (target  position  and  transmission  time  of  the 
transient).  The  rest  of  this  chapter  is  based  on  a  former  thesis  [Ref.  33]  related  to 
this  subject. 

A.  PRESENTATION  OF  TDOA  ALGORITHM 

In  this  section,  a  brief  description  of  the  TDOA  algorithm  will  be  provided  to 
the  reader  since  a  detailed  derivation  of  this  algorithm  is  not  within  the  scope  of  this 
thesis.  Before  we  continue  with  the  rest  of  our  discussion,  it  would  be  advisable  to 
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make  clear  that  the  solution  to  the  locahzation  problem,  which  is  attempted  to  be 
given  here,  is  conducted  in  two  dimensions  in  the  x  -  y  plane.  Let  us  assmne  that 
we  have  one  source  (target)  and  N  receivers  (buoys).  The  location  for  each  receiver 
is  given  by  the  vector 


ri  = 


Xi 

Vi 


(VI.l) 


where  Xi  and  yi  are  the  x-y  coordinates  of  buoy  i.  Let  the  arrival  time  of  the  transient 
at  buoy  i,  defined  as  the  time  it  takes  the  leading  edge  of  the  transient  to  travel  fi:om 
the  source  to  the  receiver,  be  denoted  as  In  a  similar  manner,  the  transmission 
time  of  the  transient  from  the  target  will  be  denoted  as  t,  and  the  desired  position  of 
the  target  is  defined  by 


r  = 


(VI.2) 


where  x  equals  the  target’s  x  coordinate  position,  y  equals  the  target’s  y  coordinate 
position. 

In  order  to  estimate  the  desired  variables,  t  and  r,  we  have  to  solve  the  following 
equation  in  matrix  form 


Ar  =  qi  +  s, 


(VI.3) 


where 


A  = 


Xi  —  X2 
X2  -  Xz 


yi-y2 

y2-y2 


XN-1  —  Xn  yN-1  —  yN 


q=  c" 


tl-t2 

t2-tz 

tN-l  —  tN 


(VI.4) 


(VI.5) 
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and 


Ilrilp-INP-Aj  +  Ai 
INP-llrslp-A^  +  AI 


(VI.6) 


The  variable  c  represents  the  sound  speed,  which  in  our  case  will  be  assumed  constant 
throughout  the  whole  computational  grid.  Equation  (VI.3)  is  simply  the  transforma¬ 
tion  of  the  nonhnear  TOA  equation 


(X(  -  xf  +  (Vi  -  yf  -  c"(ti  -tf  =  0  {VI.7) 


into  a  linear  TDOA  equation  which  relates  target  position  to  transmission  time 

2;n+l)"by(yn  J/n+l)  “  C  ^n+l)  "t"  ^  (^n  Vn+l  ^  ^n+l)’ 

(VI.8) 

This  equation  is  formed  by  subtracting  two  equations  like  (VI.7)  for  successive  pairs 
of  the  N  receivers. 

The  solution  of  Eq.  (VI.3)  is  computed  in  a  least  squares  sense  in  order  to 
find  the  target  range  r  in  terms  of  the  transmission  time  t.  Equation  (V1.3)  becomes 

r  =  gt  -4-  h,  (VI.9) 

where 

g  =  (A)+q,  (VI.10) 

h  =  (A)+s,  (VI.ll) 

and  -I-  denotes  the  pseudoinverse.  Here  we  make  the  distinction  that  when  AT  =  3 
the  expected  solution  is  exact  for  the  corresponding  set  of  equations,  but  for  the  case 
iV  >  3  the  result  is  the  least  squares  solution,  which  minimizes  the  equation  error  for 
the  same  set  of  Eqs.  (V1.7). 

After  some  calculations,  we  come  up  with  the  following  formula  for  the 
range  equation: 

|K-gt-h|P  =  c^(4-i)".  {VI.12) 
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which  after  expanding  and  simplification  produces 

at^  +  2bt  +  d  =  0, 

(VI.  13) 

where 

«=llg|^-c^ 

(VI.  14) 

6  =  •  h  -  •  Tk  +  cHk, 

(VI.  15) 

d=||h-rk|p-c^t|. 

(VI.16) 

By  solving  the  quadratic  formula  (VI. 13),  we  find  two  possible  solutions  for  transmis¬ 
sion  time  t,  but  only  one  of  them  is  correct.  We  must  use  the  correct  solution  in  Eq. 
(VI.9)  in  order  to  find  the  location  of  the  target. 

In  theory,  the  solution  for  t  is  independent  of  which  particular  range  equation 
(k)  was  used  to  develop  the  quadratic  equation.  In  practice,  because  of  measurement 
errors,  the  solutions  will  be  different.  Therefore,  it  is  advisable  to  develop  N  -  1 
estimates  for  t  using  k  =  discarding  any  erroneous  values  and  average 

the  remaining  solutions.  This  multiple  solution  technique  also  resolves  which  root  is 
correct.  The  result  can  then  be  used  in  (VI.9)  to  find  target  position. 

B.  MATLAB  IMPLEMENTATION 

This  section  explains  how  the  localization  testing  was  implemented  in  MAT- 
LAB  using  data  produced  by  the  MMPE  propagation  model.  The  problem  was  rep¬ 
resented  by  a  “target” ,  which  is  the  sound  source,  and  two  or  more  sonobuoys,  which 
are  the  receivers.  More  specifically,  we  specified  a  cartesian  coordinate  system  and 
set  the  true  position  of  the  target.  The  buoy  positions  were  given  to  the  MATLAB 
program  relative  to  the  target’s  position.  This  was  done  by  having  as  inputs  the 
quadrant  number ,  the  distance  (range)  to  the  buoy,  and  the  distance  in  x-coordinate 
to  the  buoy.  In  this  way,  it  was  possible  to  create  any  geometric  scenario  in  two 
dimensions. 
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After  specifying  the  geometry  of  the  scenario,  we  had  to  find  the  value  of 
the  sound  speed  that  would  be  used  for  each  environment.  To  be  as  realistic  as 
possible  with  each  case,  we  evaluated  the  group  speed  of  a  specific  mode  that  would 
dominate  during  propagation  of  the  signal.  This  was  executed  for  each  water  column 
corresponding  to  each  buoy  location  in  the  range-dependent  cases  and  only  once  in 
range-independent  cases,  since  the  sound  speed  profile  remains  the  same  throughout 
the  computational  grid.  For  the  range-dependent  environments,  we  averaged  those 
values  of  the  group  speed  for  each  buoy  in  order  to  end  up  with  one  value  to  use  for 
the  rest  of  the  algorithm. 

After  these  steps,  the  data  from  MMPE  (after  being  modulated/demodulated 
as  described  in  Chapter  4)  was  read  into  the  localization  program.  Then  using  all 
the  previously  described  methods  (subspace  -  classical)  we  evaluated  the  TDOAs  for 
all  combinations  of  pairs  of  buoys.  Since  we  don’t  know  absolute  TOA  for  any  of  the 
buoys,  we  set  TOA  for  the  1st  buoy  equal  to  zero,  TOA{\)  =  0,  and  measured  the 
relative  TOA  for  the  other  buoys  with  respect  to  the  first  buoy  as  a  reference. 

The  rest  of  the  algorithm  continues  as  it  was  presented  in  the  last  section  to 
evaluate  the  “Estimated  Target  Position.”  The  entire  procedure  is  repeated  for  all  the 
TDOA  methods  for  each  environment  and  for  each  case  and  using  the  same  signals. 
In  this  way,  it  is  possible  to  evaluate  the  performance  of  all  the  methods  under  the 
same  circumstances  and  compare  results  among  them.  Figmre  36  summarizes  all  the 
steps  of  the  localization  procedure. 

C.  SIMULATION  RESULTS 

In  this  section,  simulation  results  for  the  localization  problem  using  the  foirr 
MMPE  environments  are  presented.  For  each  environment,  there  are  two  cases  con¬ 
sidered.  In  the  first,  the  minimiun  nmnber  of  three  sonobuoys  is  used.  (Recall  that 
three  is  the  minimiim,  since  we  require  at  least  two  separate  TDOA  measurements  to 
locate  the  target.)  In  the  second  case,  more  than  three  sonobuoys  are  used. 
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Figure  36.  Flowchart  of  Localization  Algorithm  implemented  in  MATLAB  environ¬ 
ment  using  data  from  MMPE. 

For  each  environment  and  each  TDOA  method  we  show  a  pair  of  figures  with 
the  locations  of  the  buoys,  the  real  target  location  and  the  estimated  location  from 
each  method.  Each  figure  in  the  pair  corresponds  to  the  two  cases  mentioned  above  (3 
buoys  and  more  than  3  buoys),  and  contains  two  subfigures  with  the  results  from  the 
subspace  and  from  the  classical  methods  respectively.  Each  environment  is  accompa¬ 
nied  by  a  table.  The  tables  contain  the  error  in  target’s  position,  i.e.,  the  difference 
between  the  actual  location  and  the  estimated  one  for  each  method.  In  this  way  there 
is  a  fair  evaluation  for  all  the  methods,  since  for  this  kind  of  problem  the  difference 
of  the  two  positions  (real-estimated)  is  of  most  concern. 
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1.  Flat  Bottom 

The  characteristics  of  this  environment  are  described  in  Chapter  4.  The  true 
geometric  positions  for  two  scenarios  with  3  and  5  receivers  (buoys)  are  presented  in 
Figs.  37  and  38  respectively.  The  source  was  at  a  depth  of  150m  while  the  receivers 
were  at  a  depth  of  50m  and  ranges  of  2,  3  and  5  km  from  the  sotmce.  Table  XII  has 
the  target  position  error  for  all  methods  for  both  cases  in  this  environment. 

Observing  Table  XII,  we  conclude  that  all  methods  performed  efficiently,  but 
the  position  accuracy  was  better  for  the  subspace  methods  than  the  classical  meth¬ 
ods.  In  both  groups  of  methods  there  was  an  improvement  with  the  increase  of  the 
number  of  sonobuoys,  except  from  the  Roth  method.  The  variance  of  the  results  from 
the  subspace  methods  was  much  smaller  than  the  respective  one  from  the  classical 
methods.  Further,  the  cross-correlation  (X-Corr)  had  the  best  performance  from  the 
classical  methods  and  was  the  only  one  from  this  family  whose  results  were  closer  to 
those  of  the  subspace  methods. 


Method 

Target  Position  Error  in  (m) 

3  buoys 

5  buoys 

MUSIC 

29.68 

23.74 

Modified-MUSIC 

47.46 

24.67 

MIN-NORM 

30.16 

22.45 

PCLP 

30.16 

22.42 

ESPRIT 

30.76 

Root-MUSIC 

30.06 

21.81 

SCOT 

100.87 

67.52 

PHAT 

100.87 

67.52 

Roth 

87.17 

183.81 

X-Corr 

37.04 

31.42 

Table  XII.  Target  position  Error  for  Flat  bottom  isospeed  case. 
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Figure  38.  Localization  problem  for  Flat  bottom  isospeed  case  with  5  sonobuoys.  (a) 
Subspace  methods,  (b)  Classical  methods. 


2.  Sound  Channel 

The  characteristics  of  this  environment  are  described  again  in  Chapter  4.  The 
true  geometric  positions  of  the  two  scenarios  with  3  and  5  receivers  (buoys)  are  shown 
in  Figs.  39  and  40,  respectively.  The  somce  was  at  a  depth  of  150m  and  the  receivers 
were  at  a  depth  of  50m  and  at  ranges  of  2  and  4  km  from  the  source.  The  results  are 
also  accompanied  by  Table  XIII  presenting  the  Target  Position  Error  for  all  methods. 
Once  again  the  subspace  methods  had  better  accuracy  than  the  classical  methods 
as  a  group.  In  both  families  of  methods,  we  observe  that  increasing  the  number  of 
the  buoys  does  not,  by  itself,  assist  the  algorithm  to  more  accurately  predict  the 
target  position.  The  main  consideration  that  has  to  be  made  is  the  disposition  of 
the  receivers  and  then  the  munber  of  them.  If  the  disposition  is  suitable,  then  an 
increase  of  the  number  of  the  sonobuoys  may  improve  the  result  since,  as  we  have 
seen  in  section  A  of  this  chapter,  Eq.  VI.3  will  be  solved  either  uniquely  if  AT  =  3  or 
in  a  least  squares  sense  if  AT  >  3.  The  performance  was  similar  to  the  Flat  bottom 
case,  with  Roth  being  more  unstable  in  the  results,  and  SCOT,  PHAT  performing 
better  than  X-Corr. 


Method 

Target  Position  Error  in  (m)~| 

3  buoys 

5  buoys 

MUSIC 

19.65 

25.72 

Modified-MUSIC 

21.35 

28.42 

MIN-NORM 

19.65 

26.65 

PCLP 

19.65 

26.65 

ESPRIT 

20.50 

26.20 

Root-MUSIC 

19.65 

26.65 

SCOT 

47.04 

53.02 

PHAT 

47.04 

53.02 

Roth 

24.70 

487.13 

X-Corr 

45.92 

65.09 

Table  XIII.  Target  position  Error  for  Sound  Channel  case  . 
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Figure  39.  Localization  problem  for  Sound  Channel  case  with  3  sonobuoys.  (a) 
Subspace  methods,  (b)  Classical  methods. 
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3.  Shelf  Break 

Using  the  characteristics  for  this  environment  discussed  in  Chapter  4,  the  true 
geometric  positions  of  two  scenarios  with  3  and  6  receivers  (buoys)  are  presented  in 
Figs.  41  and  42  respectively.  The  source  was  at  a  depth  of  250m  and  the  receivers 
were  at  a  depth  of  50m  and  at  range  of  3  km  from  the  source,  located  upslope  and 
downslope  concurrently  (see  similar  Fig.  26).  Table  XIV  summarizes  the  results  of 
the  performance  for  all  methods.  In  this  range-dependent  environment,  we  see  that 
subspace  methods  have  similar  accuracy  in  their  results  compared  with  the  one  in 
the  range-independent  environments.  In  comparison  with  the  subspace  methods,  the 
classical  methods  performed  much  poorer.  This  could  be  expected  if  we  recall  their 
quality  of  results  in  the  estimation  of  TDOA  (see  Table  X.)  In  general,  there  was  an 
improvement  of  the  results  for  all  methods  (except  X-Corr)  with  an  increase  in  the 
number  of  the  sonobuoys.  This  improvement  was  more  drastic  for  Roth,  something 
that  also  was  expected  considering  the  diversity  of  the  results  from  this  method.  An¬ 
other  remark  is  that  because  of  the  complexity  of  the  range-dependent  environments 
the  buoys  disposition  is  very  critical.  An  inappropriate  disposition  can  create  larger 
errors  in  target  position  than  the  range-independent  environments. 


Method 

Target  Position  Error  in  (m) 

3  buoys 

6  buoys 

MUSIC 

24.65 

22.48 

Modified-MUSIC 

31.04 

22.64 

MIN-NORM 

25.32 

20.45 

PCLP 

25.32 

20.48 

ESPRIT 

27.29 

21.96 

Root-MUSIC 

27.29 

21.89 

SCOT 

434.35 

432.77 

PHAT 

434.35 

432.77 

Roth 

787.43 

418.26 

X-Corr 

93.53 

215.67 

Table  XIV.  Target  position  Error  for  Shelf  Break  case  . 
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Figure  41.  Localization  problem  for  Shelf  Break  case  with  3  sonobuoys.  (a)  Subspace 
methods,  (b)  Classical  methods. 
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4.  Internal  Waves 

Using  the  characteristics  for  this  environment  described  in  Chapter  4,  the  true 
geometric  positions  of  two  scenarios  are  shown  with  3  and  6  receivers  (buoys)  in 
Figs.  43  and  44,  respectively.  The  source  was  at  a  depth  of  150m  and  the  receivers 
were  at  a  depth  of  50m  and  at  range  of  2  km  from  the  source  located  on  both  the 
sides.  We  have  to  mention  that  some  of  the  receivers  were  located  in  the  section 
of  the  computational  grid  with  sinusoidal  and  soUton  perturbations  and  the  rest 
of  them  were  in  the  other  section  with  sinusoidal  perturbations  only  (see  similar 
Fig.  31).  Table  XV  contains  the  results  of  the  performance  for  all  methods.  This 
environment  is  the  most  complex  that  we  used  so  far  and  it  deviates  considerably  from 
the  basic  concept  of  the  localization  algorithm,  which  assumes  a  constant  value  for  the 
sound  speed.  For  this  environment  all  the  methods  have  decreased  performamce.  The 
difference  between  subspace  and  classical  was  the  same  as  before,  shown  in  Table  XL 
The  only  difference  is  that  the  increase  in  the  number  of  the  sonobuoys  with  different 
disposition  had  worse  results  than  in  the  Shelf  Break  environment  and  that  X-Corr 
gave  similar  results  to  the  subspace  methods. 


Method 

Target  Position  Error  in  (m)  || 

3  buoys 

6  buoys 

MUSIC 

65.90 

160.52 

Modified-MUSIC 

64.64 

153.58 

MIN-NORM 

63.02 

164.64 

PCLP 

62.29 

165.85 

ESPRIT 

63.81 

164.78 

Root-MUSIC 

66.11 

SCOT 

423.57 

849.08 

PHAT 

423.57 

849.08 

Roth 

226.81 

515.59 

X-Corr 

49.22 

160.37 

Table  XV.  Target  position  Error  for  Internal  Waves  case  . 
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Figure  44.  Localization  problem  for  Internal  Waves  case  with  6  sonobuoys.  (a) 
Subspace  methods,  (b)  Classical  methods. 


Prom  the  corresponding  figures  and  tables  for  each  environment,  one  can  con¬ 
clude  that  in  most  cases  all  the  methods  performed  quite  successfully.  The  subspace 
methods  were  more  consistent  in  their  results  as  a  group;  this  is  something  that  didn’t 
happen  with  the  classical  methods.  Further,  with  different  realizations  of  the  noise, 
the  variance  of  the  results  of  the  subspace  methods  was  much  smaller  than  for  the 
classical  methods.  The  location  accuracy  in  the  range-dependent  environments  was 
less  than  in  the  range-independent  environments  and  was  reduced  even  more  when 
the  distance  between  receivers  and  source  became  greater  (4  to  8km).  Another  major 
factor  for  the  successful  localization  of  the  target  was  the  geometrical  disposition  of 
the  receivers  relative  to  the  source  position.  In  other  words,  if  the  receivers  were 
located  in  such  a  way  that  the  bearing  from  each  pair  formed  an  angle  of  between 
60"  and  120°,  then  the  estimate  of  their  crossing  will  have  less  error  than  otherwise. 

A  larger  number  of  receivers  does  not  always  contribute  to  more  accurate 
detection.  This  can  be  justified  from  the  fact  that  when  N  is  equal  to  3,  the  solution 
from  Eq.  (VI.9)  is  unique,  while  a  value  of  AT  >  3  requires  a  least  squares  solution 
to  minimize  the  equation  error  in  the  over  specified  set  of  equations.  As  discussed  in 
Section  A  of  this  chapter,  the  solution  for  t  is  independent  of  which  particular  range 
equation  (k)  is  used  to  develop  the  quadratic  equation,  and  in  order  to  overcome 
the  obstacle  of  measurement  errors,  we  have  to  develop  iV  —  1  estimates  for  t  using 
k  =  1,2,  -  •  •  ,N  —  1,  discarding  any  erroneous  values  and  averaging  the  remaining 
solutions.  A  larger  value  of  N  does  not  always  lead  to  a  better  result.  In  other  words 
for  more  accurate  results,  the  placement  of  the  buoys  is  more  important  than  the 
number  of  the  buoys. 
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VII.  TRACKING  PROBLEM 


A.  DOPPLER  IMPLEMENTATION  IN  MMPE 

In  actual  underwater  acoustics  probleros,  the  sotuce  and/or  the  receiver  are 
not  stationary,  but  they  are  in  motion.  This  motion  causes  temporal  fluctuations, 
something  we  have  to  take  into  consideration  if  we  want  to  find  out  how  the  under¬ 
water  acoustic  channel  influences  various  apphcations  such  as  communications  [Ref. 
34].  In  om  case  we  will  consider  motion  only  due  to  the  source  (target),  and  the 
receivers  (buoys)  will  be  stationary.  Figure  45  gives  a  characteristic  example  of  the 
effect.  The  result  of  the  source  motion  is  to  provoke  a  change  in  the  value  of  the  actual 
source  frequency  observed,  fsi&),  in  other  words  in  the  received  frequency  at  buoy  i 
compared  with  the  original  transmitted  frequency,  /r-  The  received  frequency,  fs(^), 
will  depend  on  the  transmitted  frequency,  fr,  the  observation  angle,  6,  the  somce 
speed.  Vs,  and  the  direction  of  motion,  05  [Ref.  35], 


fsiO)  = 


_ fr _ 

1  -  cos(«  -  it's) 


(vn.i) 


I 


Figure  45.  Effect  of  Doppler  due  to  linear  motion  of  the  soiuce. 
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According  to  [Ref.  35]  when  |<^s|  <  90°  then  the  observed  frequency  bandwidth 
downrange,  will  be  equal  to 

BW^S^  =  BWt  +  ^  (/r,max  +  fT,rrnn\ Sm ,  (VIL2) 

where  BWt  is  the  actual  transmitted  bandwidth,  and  fx^max  and  /r.min  are  the  actual 
transmitted  maximum  and  minimum  frequencies,  respectively,  of  the  transmitted 
band.  The  observed  center  frequency  downrange  is  defined  as  [Ref.  35] 

r-i“  '^S 

/s,c  =  /r,c  +  ^  (/r,  max  fTymin  |sin0s|),  (VII.3) 

where  fT,c  is  the  actual  transmitted  center  frequency.  When  |05|  >  90°,  then  the 
observed  bandwidth  downrange  will  be  [Ref.  35] 

=  BWt  +  y  (/r.maxi  sin <^5]  +  ,  (VII.4) 

and  the  observed  center  frequency  is  equal  to 

fs7c  =  Me  +  —  iM,max\  SUl  +  M,min)  •  (VII.5) 

C/ 

The  equation  for  the  PE  starting  field  in  the  vertical  wavenumber  domain  ba-<s 
the  general  form  [Ref.  35] 


^{r  =  0,  k„  f)  =  f)  -  /),  (VII.6) 


where  the  fxmetions  'ipo{kz,  f)  represent  the  starting  field  in  free  space.  The  role  of  the 
exponentials  in  Eq.  (VII.6)  is  the  creation  of  the  appropriate  interference  structxnre 
between  the  source  (target),  which  hes  at  depth  zs  and  its  image  about  the  plane  of 
the  free  surface  (z  =  0).  The  vertical  wavenumber  is  defined  by 


kz  =  ko  sin  $  = 


(VII.7) 


where  6  is  the  angle  of  propagation  relative  to  horizontal  and  ko  is  the  reference 
wavenumber. 
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For  the  rest  of  this  chapter,  we  assume  that  we  have  a  point  source,  which  is 
modeled  in  such  way  that  [Ref.  35] 


^{kzj)  =  0!{ko)S{f), 


(VII.8) 


where 


Oi(^ko') 


I  iRo 

2'Kko' 


(VIL9) 


and  S{f)  is  the  amphtude  spectrum  of  the  pulse  in  the  frequency  domain.  The  MMPE 
uses  a  Hanning  window  m  order  to  create  the  aforementioned  spectrum,  defined  by 

cos^  ,  1/  -  /cl  <  ^ 


S{f)  = 


0, 


BW 


l/-/o|>  2 


(VII.IO) 


where  /c  is  the  center  frequency  of  the  transmitted  band  and  is  the  transmission 
bandwidth. 

Now  that  there  is  motion  of  the  source,  Eq.  (VII.6)  should  be  stated  in 
accordance  with  the  transmission  parameters  (see  [Ref.  35])  as 


i’{r  =  0,k„f)^  fy)  -  h). 


(VII.ll) 


where 


kz,T  = 


ko  sin  $ 

l  +  ^cos{d-cj)s)  ■ 


(VII.  12) 


The  frequency  in  the  enviromental  frame  is  defined  by  the  transmitted  frequency,  /r. 


with  the  assistance  of  Eq.  (VII.  1),  which  determines  the  relationship  between  /x  and 
fs{0).  Finally  the  source  spectrum  S{f)  will  take  the  form  [Ref.  35] 


Sifr)  =  -5 


/ 


i  +  ^cos{e -4>s) 


(VII.  13) 


Finally  the  factor  a  has  to  be  replaced  by  the  corresponding  transmitted  value  ay. 

One  of  the  diflFerences  between  the  two  cases,  with  and  without  Doppler,  is  the 
different  size  of  the  observed  frequency  bandwidth.  In  order  to  be  able  to  compare 
results,  the  computational  bandwidth  of  the  non-Doppler  case  has  to  be  set  to  the 
same  size  as  the  Doppler  case.  This  is  achieved  by  applying  a  weighting  of  zero  to 
those  frequency  bins  that  are  outside  of  the  actual  non-Doppler  bandwidth. 
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B.  MATLAB  IMPLEMENTATION  AND  SIMULATION 
RESULTS 

For  the  tracking  problem,  the  MATLAB  implementation  is  almost  the  same 
as  the  implementation  for  the  locsdization  problem.  However,  there  are  two  mflin 
dijfferences.  The  first  is  the  different  manipulation  of  the  MMPE  data,  since  the  code 
has  now  been  changed  in  order  to  give  us  the  effect  of  the  Doppler.  The  second  is  an 
addition  to  the  previous  code,  which  has  as  inputs  the  coordinates  of  the  buoys,  the 
estimated  coordinates  of  the  source,  the  received  signals  corresponding  to  each  buoy, 
and  the  evaluated  group  speed  (see  Chapter  6,  section  2).  The  computed  outputs  will 
then  be  the  course  and  the  speed  of  the  target.  Regarding  the  new  manipulation 
of  the  MMPE  data,  most  of  the  procedure  was  shown  in  the  previous  section.  So 
after  re-mapping  the  original  “MMPE”  signal  with  the  desired  frequency  resolution 
to  accotmt  for  the  shift  of  the  center  frequency,  we  take  the  product  of  the  real 
part  of  the  re-mapped  signal  with  the  phasor  where  fc  is  the  center 

frequency  of  the  received  signal,  and  Time  is  the  re-mapped  time  vector  with  the 
desired  resolution.  Finally  there  is  an  interpolation  of  the  pressure  vector  from  its 
original  time  spanning  vector  Time  to  a  fixed  time  spanning  vector  Timeax  which 
is  the  same  for  aU  received  signals. 

For  the  computation  of  the  course  and  the  speed  of  the  target,  the  procedure 
that  we  implemented  was  the  following.  First,  we  find  the  spectrum  of  each  received 
signal  and  from  the  spectrum  we  extract  the  center  frequency.  Afterwards,  with  the 
coordinates  of  the  buoy  locations  and  the  estimated  target  location,  we  extract  the 
angles  between  the  relative  position  vectors  for  each  pair  of  the  buoys.  Then  we  find 
the  relative  Doppler  frequency  that  corresponds  to  each  pair  of  center  frequencies  of 
the  received  signals.  Finally  we  estimate  the  course  of  the  target  for  each  buoy  by  the 
following  equations.  Note  that  for  the  simplicity  of  the  problem  we  iised  only  three 
receivers  (N=3),  otherwise  the  system  of  equations  becomes  very  complex. 

Course{\)  =  acos(A(l)) -f  B(l)  (VII.14) 
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Course(2)  =  acos{A{2))  +  B{2)  (VII.15) 

Course(2t)  =  acos{A(3))  +  B{S)  (VII.16) 


where 

A(l) 

(iq + Ki)  -  Ki 
iK^K, 

(VII.  17) 

^(2) 

{KI  +  KI)  -  Kl 

2K3K1 

(VII.  18) 

A(3) 

{Kl  +  KI)  -  KI 

2K1K2 

(VII.  19) 

and 

B(l) 

=  i  (*(1)  +  ,^(2)) 

(VII.20) 

B{2) 

=  2  ■*' 

(VII.21) 

B{3) 

=  5(^(3) +  .#.(1)) 

(VII.22) 

and 

= 

A/i2ci 

sin  (1(0(1)  -  0(2))) 

(VII.23) 

K2  = 

A/23ci 

sin  (1(0(2)  -  0(3))) 

(VII.24) 

A/3ici 

sin  (5(0(3)  -  0(1))) 

(VII.25) 

In  Eqs.  (VII.20)  -  (VII.22)  and  (VII.23)  -  (VII.25),  the  quantities  ^(1),  0(2),  0(3) 
represent  the  angles  in  the  horizontal  plane  of  the  position  vectors  of  the  buoys  and  in 
Eqs.  (VII.23)  -  (VII.25)  the  variables  A/12,  A/23,  A/31  stand  for  the  relative  Doppler 
frequencies  between  the  center  frequencies  of  the  received  signals.  In  the  end,  the 
course  of  the  target  will  be  the  average  of  the  outcomes  of  Eqs.  (V1I.14)  -  (VIL16). 

Regarding  the  speed,  it  will  be  the  average  of  the  following  calculations 


Speed{l) 


_ A/12C _ 

2/cSin  sin  (Course{l)  — 


(VII.26) 
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Speed{2) 

Speed{3) 


_ A/23C _ 

2/,  sin  sin  (Course(2)  - 

_ A/31C _ 

2/cSin(^/W)  sin  (Course{3)  — 


(VII.27) 

(VIL28) 


where  fc  is  the  theoretical  center  frequency  of  the  transmitted  signal.  Since  this 
quantity  is  unlcnown  in  our  case,  we  substitute  it  with  the  average  of  the  center 
frequencies  of  the  received  signals.  This  substitution  will  not  significantly  distort 
the  result,  since  the  shift  of  the  center  frequency  due  to  Doppler  is  not  large  enough 
(low  speeds)  and  the  difference  between  the  fc  (transmitted)  and  the  average  of  fc 
(received)  is  small.  Figure  46  gives  a  geometrical  overview  of  the  problem. 


Figmre  46.  Top  view  of  the  Tracking  problem. 

Some  simulation  results  will  follow,  one  from  each  environment.  For  each 
case  there  will  be  figures  indicating  the  geometric  scenario  and  the  solution  that  each 
method  gives  for  the  locahzation  problem  and  also  tables  with  Target’s  Position  Error 
in  (m)  and  the  Target’s  motion  Data,  Course{0  -  359°)  and  Speed{m/s). 

1.  Flat  Bottom 

The  characteristics  of  this  range-independent  environment  remain  the  same  as 
were  presented  in  Chapters  4,  5  and  6.  In  this  case,  the  source  is  at  depth  100m  and 
the  3  buoys  are  located  at  ranges  2,  3  and  5km  and  at  depth  30m  as  shown  in  Fig. 
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47.  The  source  course  and  speed  are  40"  and  15m/s,  respectively.  Table  XVI  presents 
the  results  for  the  localization  and  tra.cking  accuracy. 


Method 

Target’s  Position  Error  in  (m) 

MUSIC 

130.8 

Modified-MUSIC 

120.2 

MIN-NORM 

132.0 

PCLP 

131.1 

ESPRIT 

137.4 

Root-MUSIC 

137.1 

SCOT 

1010.0 

PHAT 

1010.0 

Roth 

1293.0 

X-Corr 

53.6 

Target’s  Motion 

Estimated 

Actual 

Course{0  —  359°) 

37.73 

40.00 

Speed  (m/s) 

13.95 

15.00 

Table  XVI.  Target’s  Motion  h  position  Error  for  Flat  bottom  isospeed  case  (Number 
of  receivers  =  3). 
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2.  Sound  Channel 

As  before,  this  range-independent  environment  retains  the  same  properties  as 
shown  in  Chapters  4,  5  and  6.  For  this  case,  the  source  is  at  depth  150m  and  the  3 
buoys  are  located  at  ranges  2  and  4km  and  at  depth  40m  as  shown  in  Fig.  48.  The 
source  course  and  speed  are  75“  and  lOm/s,  respectively.  Table  XVII  presents  the 
results  for  the  localization  and  tracking  accuracy. 


Method 

Target’s  Position  Error  in  (m) 

MUSIC 

126.0 

Modified-MUSIC 

131.7 

MIN-NORM 

127.6 

PCLP 

126.3 

ESPRIT 

130.6 

Root-MUSIC 

129.4 

SCOT 

2234.5 

PEAT 

2234.5 

Roth 

4413.0 

X-Corr 

229.3 

Target’s  Motion 

Estimated 

Actual 

Course  {0  -  3'59“) 

75.17 

75.00 

Speed  (m/s) 

11.41 

10.00 

Table  XVII.  Tsurget’s  Motion  &  position  Error  for  Sovmd  Channel  case  (Number  of 
receivers  =  3). 
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3.  Shelf  Break 

The  properties  of  the  current  range-dependent  environment  are  the  same  as 
described  in  previous  Chapters.  In  this  case,  the  source  is  at  depth  250m  and  the  3 
buoys  are  located  at  ranges  3  and  4km  and  at  depth  100m  as  shown  in  Fig.  49.  The 
source  course  and  speed  are  100®  and  8m/s,  respectively.  Table  XVIII  presents  the 
results  for  the  localization  and  tracking  accuracy. 


Method 

Target’s  Position  Error  in  (m) 

MUSIC 

259.2 

Modified-MUSIC 

289.3 

MIN-NORM 

282.8 

PCLP 

275.6 

ESPRIT 

335.7 

Root-MUSIC 

303.5 

SCOT 

671.4 

PHAT 

671.4 

Roth 

1947.6 

X-Corr 

185.2 

Teirget’s  Motion 

Estimated 

Actual 

Course  (0  -  359®) 

105.64 

100.00 

Speed  (m/s) 

8.43 

8.00 

Table  XVIII.  Target’s  Motion  &  position  Error  for  Shelf  Break  case  (Number  of 
receivers  =  3). 
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Figure  49.  Tracking  problem  for  Shelf  Break  case  with  N=3  sonobuoys.  (a)  Subspace 


methods,  (b)  Classical  methods. 
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4.  Internal  Waves 

Finally  this  range-dependent  environment  retains  the  properties  that  we  have 
already  seen  in  our  previous  discussion.  During  the  present  simulation,  the  somce  is 
at  depth  150m  and  the  positions  for  the  3  buoys  are  at  range  2km  and  at  depth  50m 
as  shown  in  Fig.  50.  The  somce  comse  and  speed  are  230"  and  5m/s,  respectively. 
Table  XIX  presents  the  results  for  the  localization  and  tracking  accmacy. 


Method 

Target’s  Position  Error  in  (m) 

MUSIC 

176.5 

Modified-MUSIC 

186.3 

MIN-NORM 

176.7 

PCLP 

176.2 

ESPRIT 

181.4 

Root-MUSIC 

181.3 

SCOT 

938.5 

PHAT 

938.5 

Roth 

1504.8 

X-Corr 

103.2 

Target’s  Motion 

Estimated 

Actual 

Course{0  —  359") 

228.77 

230.00 

Speed  (m/s) 

4.21 

5.00 

Table  XIX.  Target’s  Motion  &  position  Error  for  Internal  Waves  case  (Number  of 
receivers  =  3). 
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Figure  50.  Tracking  problem  for  Internal  Waves  case  with  N=3  sonobuoys.  (a) 
Subspace  methods,  (b)  Classical  methods. 


Prom  the  presented  results,  we  observe  that  the  accuracy  of  the  localization  of 
the  target  has  been  decreased,  compared  with  the  cases  where  there  was  no  movement 
of  the  source  (zero  Doppler).  Especially  for  the  classical  methods,  the  error  in  posi¬ 
tion’s  target  was  large  enough  to  say  that  they  failed  to  estimate  the  location  of  the 
source.  One  possible  reason  for  the  poor  performance  of  the  methods  on  the  target 
localization  is  the  approximation  that  we  made  during  the  manipulation  of  the  data 
from  MMPE  when  we  were  trying  to  set  all  signals  in  a  fixed  time  spanning  vector. 
On  the  other  hand,  the  estimation  of  the  target  course  and  speed  was  quite  accurate, 
or  at  least  inside  the  acceptable  limits  of  error,  ±5®  for  the  course  and  ±2m/s  for  the 
speed. 
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VIII. 


CONCLUSIONS  AND  FUTURE  WORK 


A.  CONCLUSIONS 

In  this  thesis,  new  techniques  were  developed  to  estimate  the  Time  Difference 
of  Arrival  (TDOA)  between  two  received  signals  at  two  separate  locations  which 
contain  the  same  transient  but  different  noise.  These  techniques  are  based  on  subspace 
methods,  a  class  of  techniques  based  on  the  concept  of  signal  and  noise  subspaces 
associated  with  the  correlation  matrix  for  a  random  process,  and  whose  principal 
application  is  to  locate  narrowband  hues  in  a  spectrum  or  to  estimate  bearing  of 
narrowband  sources  using  array  processing.  As  far  as  we  know,  the  use  of  subspace 
methods  to  estimate  TDOA  for  broadband  transient  sources  is  new.  Our  objectives 
were  to  implement  the  subspace  methods  in  our  problem,  and  to  test  them  thoroughly 
comparing  their  performance  with  traditional  methods  of  TDOA  based  on  generalized 
cross-correlation  (“classical  methods”).  The  testing  was  conducted  using  synthetic 
data  generated  in  MATLAB  and  data  from  an  acoustic  propagation  model  (MMPE). 
Tests  were  conducted  under  different  scenarios  using  various  transients,  noise,  or 
environmental  conditions. 

The  thesis  went  further  to  apply  these  methods  in  the  localization  of  a  target 
using  the  TDOA  ranging  algorithm  and  data  from  MMPE  generated  for  various  ocean 
environments.  In  this  way,  we  tried  to  simulate  realistic  conditions  in  the  ocean  as 
much  as  possible.  This  second  step  was  further  expanded  by  the  use  of  Doppler  data 
from  the  MMPE  algorithm  to  estimate  not  only  the  position  but  also  course  and 
speed;  in  other  words,  to  track  the  source. 

Prom  the  various  types  of  testing,  we  can  conclude  that  the  subspace  methods 
provide  good  results  in  almost  all  cases.  As  a  group,  the  subspace  methods  were 
more  consistent  in  the  accuracy  of  the  results  than  the  classical  methods.  This  was 
observed  not  only  in  the  TDOA  testing  but  also  in  the  locahzation  and  tracking 
experiments,  where  the  subspace  methods  were  more  consistently  accurate  than  the 
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classical  methods. 

The  plots  produced  by  the  subspace  methods  provide  a  clear  and  more  easily 
identified  correct  peak  than  those  corresponding  to  the  classical  methods  and  were 
not  mfluenced  as  much  by  the  noise.  In  addition  some  of  these  subspace  methods 
can  also  be  used  to  produce  direct  numerical  estimates  for  the  time  delay  without  the 
need  to  search  for  a  peak  (ESPRIT  -  root  MUSIC). 

Almost  aU  methods  (classical  and  subspace)  performed  satisfactorily  in  the 
localization  problem.  However  accuracy  tended  to  be  more  degraded  for  the  range- 
dependent  problems.  This  is  not  too  surprising  since  the  basic  localization  algorithm 
makes  the  assumption  of  a  constant  speed  of  sound.  In  the  tracking  experiments, 
the  subspace  methods  had  hmited  performance  for  estimation  of  target  position  but 
much  better  performance  for  estimation  of  comse  and  speed. 

B.  SUGGESTIONS  FOR  FUTURE  DEVELOPMENT 

The  results  of  this  thesis  show  promising  results  for  the  subspace  methods, 
and  their  apphcation  to  target  localization  and  tracking  via  TDOA  estimation.  The 
testing  should  be  continued  with  more  complex  environments  using  acoustic  models 
and  also  with  real  data  from  ocean  experiments.  A  comparison  of  the  subspace 
methods  with  other  families  of  methods  for  TDOA  (e.g.,  bicorrelation  and  bispectra 
[Ref.  6,  7,  8,  9]),  or  even  a  combination  of  these,  would  give  a  wider  view  of  the 
possible  solutions  to  the  problem.  Finally,  implementation  of  a  three-dimensional 
locahzation  algorithm  may  give  more  realistic  results  and  could  assist  in  the  tracking 
problem,  provided  that  it  can  overcome  obstacles  such  as  multipath  propagation. 


APPENDIX  A.  SYNTHETIC  DATA 


This  appendix  presents  sample  figures  and  tables  with  summary  results  of  the 
time  delay  estimation  obtained  from  various  forms  of  synthetic  data.  The  transients 
that  were  used  are  the  following: 

•  Exponential 

•  Sinusoidal 

•  Modified-Ejqjonential 

•  Chirp 

•  Modified-Chirp 

The  definitions  of  the  transients  are  given  in  Chapter  4  Eqs.  (IV.  2)  -  (IV.  8). 
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Method 


MUSIC 


Modified  MUSIC 


MIN-NORM 


ESPRIT 


Root-MUSIC 


SCOT 


PHAT 


Roth 


X-Corr 


Parameters 


Cov-Mat:  10 


Cov-Mat:  20 


Cov-Mat:  10 


Cov-Mat:  20 


Cov-Mat:  10 
Cov-Mat:  20 


Cov-Mat:  10 


Cov-Mat:  10 


Cov-Mat:  20 


Cov-Mat:  10 


Cov-Mat:  20 


#-Segs:  1 


#-Segs:  2 


#-Segs:  1 
#-Segs:  2 


#-Segs:  1 
#-Segs:  2 


#-Segs:  1 


#-Segs:  2 


SNR  15  dB 

SNR  10  dB 

SNR  5  dB 

101 

100 

106 

99 

99 

104 

101 

100 

107 

99 

99 

104 

100 

100 

no 

99 

98 

102 

100 

100 

109 

99 

98 

103 

100 

101 

107 

99 

98 

103 

101 

100 

106 

99 

99 

104 

102 

385 

5 

101 

102 

9 

102 

385 

5 

101 

202 

98 

-228 

315 

233 

101 

201 

9 

98 

107 

93 

101 

102 

109 

Table  XX.  Exponential  transient:  length  L=20s;  white  noise  =  2):  actual 
TDOA=100s. 
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Figrire  51.  Exponential  transient:  Actual  Tdoa=100s;  SNR=15dB;  White-Noise  (cTp). 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg- 
ments=l. 
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(b) 

Figure  52.  Exponential  transient:  Actual  Tdoa=100s;  SNR=15dB;  White-Noise  (a^) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 
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Figure  53.  Exponential  transient:  Actual  Tdoa=100s;  SNR=10dB;  White-Noise  (a^) 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg 
ments=l. 


Figure  54.  Exponential  transient:  Actual  Tdoa=100s;  SNR=10dB;  White-Noise  (a^) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 


Figure  55.  Exponential  transient:  Actual  Tdoa=100s;  SNR=05dB;  White-Noise  (cr^) 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg 
ments=l. 


0 


MinNorm:  102 


X-Corr  109 


Figure  56.  Exponential  transient:  Actual  Tdoa=100s;  SNR=05dB;  White-Noise  (a^) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 


Method 

Parameters 

SNR  15  dB 

SNR  10  dB 

SNR  5  dB 

MUSIC 

Cov-Mat:  10 

-10 

-10 

-10 

-8 

-10 

Modified  MUSIC 

Cov-Mat:  10 

-10 

-10 

-10 

Cov-Mat:  20 

-11 

-8 

-10 

MIN-NORM 

Cov-Mat:  10 

-10 

-10 

-9 

Cov-Mat:  20 

-10 

-9 

-10 

PCLP 

Cov-Mat:  10 

-10 

-10 

-9 

Cov-Mat:  20 

-10 

-9 

-10 

ESPRIT 

Cov-Mat:  10 

-10 

-10 

-9 

Cov-Mat:  20 

-10 

-9 

-14 

Root-MUSIC 

Cov-Mat:  10 

-10 

-10 

-9 

Cov-Mat:  20 

-10 

-8 

-10 

SCOT 

#-Segs:  1 

-10 

-10 

-76 

#-Segs:  2 

-10 

-11 

PHAT 

-10 

-10 

-76 

#-Segs:  2 

-10 

-15 

-91 

Roth 

#-Segs:  1 

-371 

-42 

-131 

#-Segs:  2 

-10 

-209 

169 

X-Corr 

#-Segs:  1 

-10 

-9 

-9 

#-Segs:  2 

-10 

-11 

-5 

Table  XXI.  Exponential  transient:  length  L=50s;  white  noise  (a^  =  2);  actual 
TDOA=-10s. 
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■•■(sec)  T(seo) 

(b) 

Figure  57.  Exponential  transient:  Actual  Tdoa=-10s;  SNR=15dB;  White-Noise  (cr^) 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg 
ments=l. 


Figure  58.  Exponential  transient:  Actual  Tdoa=-10s;  SNR=15dB;  White-Noise  ((j^) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 
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Figure  59.  Exponential  transient:  Actual  Tdoa=-10s;  SNR=10dB;  White-Noise  (<7^) 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg 
ments=l. 


Figure  60.  Exponential  transient:  Actual  Tdoa=-10s;  SNR=10dB;  White-Noise  (cr^) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 
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Figure  61.  Exponential  transient;  Actual  Tdoa^ 
(a)  Subspace  methods,  Covariance  Size=10.  (1 
ments=l. 


-10s;  SNR=05dB;  White-Noise  (a^) 
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Figure  62.  Exponential  transient:  Actual  Tdoa=-10s;  SNR=05dB;  White-Noise  (al) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg- 
ments=2. 


Method 

Pareimeters 

Damped 

Sin 

Chirp 

Damped 

Chirp 

MUSIC 

Cov-Mat:  10 

50 

49 

50 

50 

50 

Cov-Mat:  20 

50 

51 

50 

50 

50 

Modified-MUSIC 

Cov-Mat:  10 

50 

50 

50 

50 

50 

Cov-Mat:  20 

50 

50 

50 

50 

50 

MIN-NORM 

Cov-Mat:  10 

50 

49 

50 

50 

50 

Cov-Mat:  20 

50 

51 

50 

PCLP 

Cov-Mat:  10 

50 

49 

50 

50 

50 

Cov-Mat:  20 

50 

51 

50 

ESPRIT 

Cov-Mat:  10 

49 

49 

50 

50 

Cov-Mat:  20 

50 

50 

Root-MUSIC 

Cov-Mat:  10 

50 

49 

50 

50 

50 

Cov-Mat:  20 

50 

51 

50 

50 

n  50 

SCOT 

#-Segs:  1 

50 

50 

50 

50 

50 

#-Segs:  2 

50 

-74 

50 

50 

PHAT 

#-Segs:  1 

50 

50 

50 

50 

50 

#-Segs:  2 

50 

-74 

50 

50 

50 

Roth 

#-Segs:  1 

-299 

-61 

50 

50 

50 

#-Segs:  2 

50 

50 

50 

50 

X-Corr 

#-Segs:  1 

50 

50 

50 

50 

50 

#-Segs:  2 

50 

48 

50 

50 

50 

Table  XXII.  Tiransients:  same  length;  white  noise  =  2);  SNR=15dB;  actual 
TDOA=50s. 
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Figure  63.  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=15dB;  White-Noise  (al) 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg 
ments=l. 


Figure  64.  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=15dB;  White-Noise  (a^) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 


Figiire  65.  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=10dB;  White-Noise  (cr^) 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg 
ments=l. 


Figure  66.  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=10dB;  White-Noise  {al) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 
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Figure  67.  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=05dB;  White-Noise  {a^) 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg 
ments=l. 


Figure  68.  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=05dB;  White-Noise  (a^) 
(a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg 
ments=2. 


Method 

Pzirameters 

Exp 

Chirp 

Damped 

Chirp 

MUSIC 

Cov-Mat:  10 

51 

51 

50 

53 

49 

Cov-Mat:  20 

49 

51 

50 

50 

50 

Modified-MUSIC 

Cov-Mat:  10 

51 

52 

50 

53 

49 

51 

50 

50 

50 

MIN-NORM 

52 

50 

54 

49 

Cov-Mat:  20 

49 

51 

50 

50 

50 

PCLP 

■tiaaaBiaBi 

52 

50 

54 

49 

51 

50 

50 

50 

ESPRIT 

Cov-Mat:  10 

51 

52 

50 

54 

49 

51 

50 

50 

50 

Root-MUSIC 

52 

50 

54 

49 

Cov-Mat:  20 

49 

51 

50 

50 

50 

SCOT 

#-Segs:  1  , 

50 

50 

50 

50 

50 

#-Segs:  2 

52 

-255 

50 

50 

PH  AT 

50 

50 

50 

50 

50 

#-Segs:  2 

52 

48 

50 

50 

50 

Roth 

#-Segs:  1 

-502 

50 

108 

50 

50 

52 

-255 

50 

50 

50 

X~Corr 

#-Segs:  1 

50 

50 

50 

50 

50 

#-Segs:  2 

50 

48 

50 

50 

50 

Table  XXIII.  Transients:  same  length;  white  noise  =  2);  SNR=10dB;  actual 
TDOA=50s. 
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Figure  69.  Damped  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=15dB;  White- 
Noise  (<7o).  (a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  num¬ 
ber  of  segments=l. 


Figure  70.  Damped  Sinusoidal  transient;  Actual  Tdoa=50s;  SNR=15dB;  White- 
Noise  (cTo).  (a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  num¬ 
ber  of  segments=2. 


Figure  71.  Damped  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=10dB;  White- 
Noise  (a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  num¬ 
ber  of  segments=l. 


Figure  72.  Damped  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=10dB;  White- 
Noise  (cTo).  (a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  num¬ 
ber  of  segments=2. 


Figure  73.  Damped  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=05dB;  White- 
Noise  (al).  (a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  num¬ 
ber  of  segments=l. 
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Figure  74.  Damped  Sinusoidal  transient:  Actual  Tdoa=50s;  SNR=05dB;  Widte- 
Noise  (cr^).  (a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  num¬ 
ber  of  segments=2. 
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Method 

Parameters 

Exp 

Sin 

Damped 

Sin 

Chirp 

MUSIC 

Cov-Mat:  10 

51 

50 

49 

53 

51 

49 

52 

51 

48 

51 

Modified-MUSIC 

Cov-Mat:  10 

52 

49 

53 

51 

Cov-Mat:  20 

50 

51 

51 

48 

51 

MIN-NORM 

Cov-Mat:  10 

50 

50 

49 

57 

50 

Cov-Mat:  20 

49 

52 

51 

50 

50 

PCLP 

Cov-Mat:  10 

52 

50 

49 

56 

50 

49 

52 

51 

50 

50 

ESPRIT 

Cov-Mat:  10 

50 

50 

49 

53 

49 

Cov-Mat:  20 

50 

52 

51 

51 

Root-MUSIC 

Cov-Mat:  10 

49 

50 

49 

52 

50 

Cov-Mat:  20 

49 

52 

51 

r  49 

51 

SCOT 

#-Segs:  1 

92 

40 

52 

50 

50 

#-Segs:  2 

47 

48 

50 

49 

PHAT 

#-Segs:  1 

92 

40 

52 

50 

50 

#-Segs:  2 

47 

46 

56 

50 

49 

Roth 

#-Segs:  1 

248 

478 

-204 

50 

50 

#-Segs:  2 

-216 

50 

49 

X-Corr 

#-Segs:  1 

53 

48 

46 

50 

50 

#-Segs:  2 

IB 

48 

50 

50 _ 1 

Table  XXIV.  Ibansients:  same  length;  white  noise  (a^  =  2);  SNR=05dB;  actual 
TDOA=50s. 
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Figure  75.  Chirp  transient:  Actual  Tdoa=50s;  SNR=15dB;  White-Noise  (cr^).  (a) 
Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg- 
ments=l. 


T(seo)  T(seo) 

(b) 

Figure  76.  Chirp  transient:  Actual  Tdoa=50s;  SNR=15dB;  White-Noise  (cr^).  (a) 
Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg- 
ments=2. 


Figure  77.  Chirp  transient:  Actual  Tdoa=50s;  SNR=10dB;  White-Noise  (cg).  (a) 
Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg- 
ments=l. 
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Figure  78.  Chirp  transient:  Actual  Tdoa=50s;  SNR=10dB;  White-Noise  (a^).  (a) 
Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg- 
ments=2. 
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Figure  79.  Chirp  transient:  Actual  Tdoa=50s;  SNR=05dB;  White-Noise  (cr^).  (a) 
Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of  seg- 
ments=l. 
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Figure  80.  Chirp  transient:  Actual  Tdoa=50s;  SNR=05dB;  White-Noise  (a^).  (a) 
Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of  seg- 
ments=2. 


Figure  81.  Damped  Chirp  transient:  Actual  Tdoa=50s;  SNR=15dB;  White-Noise 
(a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of 
segments=l. 
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Figure  82.  Damped  Chirp  transient:  Actual  Tdoa=50s;  SNR=15dB;  White-Noise 
(^o)-  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of 
segments=2. 
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Figure  83.  Damped  Chirp  transient:  Actual  Tdoa=50s;  SNR=10dB;  White-Noise 
(cr^).  (a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of 
segments=l. 
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Figure  84.  Damped  Chirp  transient:  Actual  Tdoa=50s;  SNR=10dB;  White-Noise 
(al).  (a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of 
segments=2. 
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Figure  85.  Damped  Chirp  transient:  Actual  Tdoa=50s;  SNR=05dB;  White-Noise 
(al).  (a)  Subspace  methods,  Covariance  Size=10.  (b)  Classical  methods,  number  of 
segments=l. 
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Figure  86.  Damped  Chirp  transient:  Actual  Tdoa=50s;  SNR=05dB;  White-Noise 
{aj).  (a)  Subspace  methods,  Covariance  Size=20.  (b)  Classical  methods,  number  of 
segments=2. 
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